/[MITgcm]/MITgcm_contrib/cg2d_bench/cg2d.F
ViewVC logotype

Diff of /MITgcm_contrib/cg2d_bench/cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by ce107, Fri May 12 21:58:05 2006 UTC revision 1.3 by ce107, Wed May 31 16:27:25 2006 UTC
# Line 1  Line 1 
1    C       $Id$    
2        SUBROUTINE CG2D        SUBROUTINE CG2D
3  C     /==========================================================\  C     /==========================================================\
4  C     | SUBROUTINE CG2D                                          |  C     | SUBROUTINE CG2D                                          |
# Line 28  C     === Global data === Line 29  C     === Global data ===
29  #include "EEPARAMS.h"  #include "EEPARAMS.h"
30  #include "PARAMS.h"  #include "PARAMS.h"
31  #include "CG2D.h"  #include "CG2D.h"
32    #if defined(USE_PAPI_FLOPS) || defined(USE_PAPI_FLIPS)
33    #if defined(PAPI_PER_ITERATION) || defined(PAPI_PER_TIMESTEP)
34    #include "PAPI.h"
35    #endif
36    #endif
37    
38  C     === Routine arguments ===  C     === Routine arguments ===
39  C     myThid - Thread on which I am working.  C     myThid - Thread on which I am working.
# Line 49  C                   or they converge at Line 55  C                   or they converge at
55  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
56  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
57        INTEGER actualIts        INTEGER actualIts
       REAL    actualResidual  
58        INTEGER bi, bj                      INTEGER bi, bj              
59        INTEGER I, J, N        INTEGER I, J, N
60        REAL    err  #ifdef USE_MIXED_PRECISION
61        REAL    errSum        REAL*8    actualResidual
62        REAL    etaN        REAL*8    err
63        REAL    etaNM1        REAL*8    errSum
64        REAL    etaNSum        REAL*8    etaN
65        REAL    beta        REAL*8    etaNM1
66        REAL    alpha        REAL*8    etaNSum
67        REAL    alphaSum        REAL*8    beta
68        REAL    sumRHS        REAL*8    alpha
69        REAL    temp        REAL*8    alphaSum
70          REAL*8    sumRHS
71          REAL*8    temp
72    #else
73          Real    actualResidual
74          Real    err
75          Real    errSum
76          Real    etaN
77          Real    etaNM1
78          Real    etaNSum
79          Real    beta
80          Real    alpha
81          Real    alphaSum
82          Real    sumRHS
83          Real    temp
84    #endif
85    
86  C--   Initialise inverter  C--   Initialise inverter
87        etaNM1              = 1. D0        etaNM1              = 1. _d 0
88    
89  C--   Initial residual calculation  C--   Initial residual calculation
90        err    = 0. _d 0        err    = 0. _d 0
91        sumRHS = 0. _d 0        sumRHS = 0. _d 0
92        DO J=1,sNy        DO J=1,sNy
93         DO I=1,sNx         DO I=1,sNx
94          cg2d_s(I,J) = 0.          cg2d_s(I,J) = 0. _d 0
95           cg2d_r(I,J) = cg2d_b(I,J) -           cg2d_r(I,J) = cg2d_b(I,J) -
96       &   ( aW2d(I  ,J  )*cg2d_x(I-1,J  )+aW2d(I+1,J  )*cg2d_x(I+1,J  )       &   ( aW2d(I  ,J  )*cg2d_x(I-1,J  )+aW2d(I+1,J  )*cg2d_x(I+1,J  )
97       &    +aS2d(I  ,J  )*cg2d_x(I  ,J-1)+aS2d(I  ,J+1)*cg2d_x(I  ,J+1)       &    +aS2d(I  ,J  )*cg2d_x(I  ,J-1)+aS2d(I  ,J+1)*cg2d_x(I  ,J+1)
# Line 92  C--   Initial residual calculation Line 112  C--   Initial residual calculation
112        actualIts      = 0        actualIts      = 0
113        actualResidual = SQRT(err)        actualResidual = SQRT(err)
114        WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual        WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual
115        IF ( actualResidual .EQ. 0. ) STOP 'ABNORMAL END: RESIDUAL 0'        IF ( actualResidual .EQ. 0. _d 0) STOP 'ABNORMAL END: RESIDUAL 0'
116    
117  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
118        DO 10 N=1, cg2dMaxIters        DO 10 N=1, cg2dMaxIters
# Line 126  C--    processes. Line 146  C--    processes.
146    
147  C--    Ten extra exchanges  C--    Ten extra exchanges
148  #ifdef TEN_EXTRA_EXCHS  #ifdef TEN_EXTRA_EXCHS
149         CALL EXCH_XY_R8( cg2d_s )         DO J=1,10
150         CALL EXCH_XY_R8( cg2d_s )            CALL EXCH_XY_R8( cg2d_s )
151         CALL EXCH_XY_R8( cg2d_s )         ENDDO
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
        CALL EXCH_XY_R8( cg2d_s )  
152  #endif  #endif
153    
154  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
# Line 155  C==    q = A.s Line 168  C==    q = A.s
168    
169  #ifdef HUNDRED_EXTRA_SUMS  #ifdef HUNDRED_EXTRA_SUMS
170  C--    Hundred extra global sums  C--    Hundred extra global sums
171         CALL GSUM_R8( temp, alpha )         DO J=1,100
172         CALL GSUM_R8( temp, alpha )            CALL GSUM_R8( temp, alpha )
173         CALL GSUM_R8( temp, alpha )         ENDDO
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
        CALL GSUM_R8( temp, alpha )  
174  #endif  #endif
175    
176         alpha = temp         alpha = temp
# Line 277  C      Now compute "interior" points. Line 193  C      Now compute "interior" points.
193         err = SQRT(err)         err = SQRT(err)
194         actualIts      = N         actualIts      = N
195         actualResidual = err         actualResidual = err
196  CcnhDebugStarts  #ifdef RESIDUAL_PER_ITERATION
197  C      WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual         WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual
198  CcnhDebugEnds  #endif
199         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTargetResidual ) GOTO 11
200         CALL EXCH_XY_R8(cg2d_r )         CALL EXCH_XY_R8(cg2d_r )
201    #ifdef PAPI_PER_ITERATION
202    #ifdef USE_PAPI_FLOPS
203           call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
204           WRITE(6,'(F10.3,A7,F10.3,A37,I8)')
205         $      mflops, ' user ', mflops*proc_time/real_time,
206         $      ' wallclock Mflop/s during iteration ', N
207    #else
208    #ifdef USE_PAPI_FLIPS
209           call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
210           WRITE(6,'(F10.3,A7,F10.3,A37,I8)')
211         $      mflops, ' user ', mflops*proc_time/real_time,
212         $      ' wallclock Mflip/s during iteration ', N
213    #endif
214    #endif
215    #endif
216     10 CONTINUE     10 CONTINUE
217     11 CONTINUE     11 CONTINUE
218        CALL EXCH_XY_R8(cg2d_x )        CALL EXCH_XY_R8(cg2d_x )
219    #ifdef PAPI_PER_TIMESTEP
220    #ifdef USE_PAPI_FLOPS
221           call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
222           WRITE(6,'(F10.3,A7,F10.3,A37,I8)')
223         $      mflops, ' user ', mflops*proc_time/real_time,
224         $      ' wallclock Mflop/s during iteration ', N
225    #else
226    #ifdef USE_PAPI_FLIPS
227           call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
228           WRITE(6,'(F10.3,A7,F10.3,A37,I8)')
229         $      mflops, ' user ', mflops*proc_time/real_time,
230         $      ' wallclock Mflip/s during iteration ', N
231    #endif
232    #endif
233    #endif
234        WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual        WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual
235    
236  C     Calc Ax to check result  C     Calc Ax to check result
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22