| 1 |  | C       $Id$ | 
| 2 | SUBROUTINE CG2D | SUBROUTINE CG2D | 
| 3 | C     /==========================================================\ | C     /==========================================================\ | 
| 4 | C     | SUBROUTINE CG2D                                          | | C     | SUBROUTINE CG2D                                          | | 
| 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. | 
| 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) | 
| 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 | 
| 297 | err = SQRT(err) | err = SQRT(err) | 
| 298 | actualIts      = N | actualIts      = N | 
| 299 | actualResidual = err | actualResidual = err | 
| 300 | CcnhDebugStarts | #ifdef RESIDUAL_PER_ITERATION | 
| 301 | C      WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual | WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual | 
| 302 | CcnhDebugEnds | #endif | 
| 303 | IF ( err .LT. cg2dTargetResidual ) GOTO 11 | IF ( err .LT. cg2dTargetResidual ) GOTO 11 | 
| 304 | CALL EXCH_XY_R8(cg2d_r ) | CALL EXCH_XY_R8(cg2d_r ) | 
| 305 |  | #ifdef PAPI_PER_ITERATION | 
| 306 |  | #ifdef USE_PAPI_FLOPS | 
| 307 |  | call PAPIF_flops(real_time, proc_time, flpops, mflops, check) | 
| 308 |  | WRITE(6,'(F10.3,A7,F10.3,A37,I8)') | 
| 309 |  | $      mflops, ' user ', mflops*proc_time/real_time, | 
| 310 |  | $      ' wallclock Mflop/s during iteration ', N | 
| 311 |  | #else | 
| 312 |  | #ifdef USE_PAPI_FLIPS | 
| 313 |  | call PAPIF_flips(real_time, proc_time, flpops, mflops, check) | 
| 314 |  | WRITE(6,'(F10.3,A7,F10.3,A37,I8)') | 
| 315 |  | $      mflops, ' user ', mflops*proc_time/real_time, | 
| 316 |  | $      ' wallclock Mflip/s during iteration ', N | 
| 317 |  | #endif | 
| 318 |  | #endif | 
| 319 |  | #endif | 
| 320 | 10 CONTINUE | 10 CONTINUE | 
| 321 | 11 CONTINUE | 11 CONTINUE | 
| 322 | CALL EXCH_XY_R8(cg2d_x ) | CALL EXCH_XY_R8(cg2d_x ) | 
| 323 |  | #ifdef PAPI_PER_TIMESTEP | 
| 324 |  | #ifdef USE_PAPI_FLOPS | 
| 325 |  | call PAPIF_flops(real_time, proc_time, flpops, mflops, check) | 
| 326 |  | WRITE(6,'(F10.3,A7,F10.3,A37,I8)') | 
| 327 |  | $      mflops, ' user ', mflops*proc_time/real_time, | 
| 328 |  | $      ' wallclock Mflop/s during iteration ', N | 
| 329 |  | #else | 
| 330 |  | #ifdef USE_PAPI_FLIPS | 
| 331 |  | call PAPIF_flips(real_time, proc_time, flpops, mflops, check) | 
| 332 |  | WRITE(6,'(F10.3,A7,F10.3,A37,I8)') | 
| 333 |  | $      mflops, ' user ', mflops*proc_time/real_time, | 
| 334 |  | $      ' wallclock Mflip/s during iteration ', N | 
| 335 |  | #endif | 
| 336 |  | #endif | 
| 337 |  | #endif | 
| 338 | WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual | WRITE(6,*) ' CG2D iters, err = ', actualIts, actualResidual | 
| 339 |  |  | 
| 340 | C     Calc Ax to check result | C     Calc Ax to check result |