/[MITgcm]/MITgcm/model/src/cg3d.F
ViewVC logotype

Diff of /MITgcm/model/src/cg3d.F

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

revision 1.24 by jmc, Wed Jun 8 01:46:34 2011 UTC revision 1.25 by jmc, Fri May 11 23:34:06 2012 UTC
# Line 14  CBOP Line 14  CBOP
14  C     !ROUTINE: CG3D  C     !ROUTINE: CG3D
15  C     !INTERFACE:  C     !INTERFACE:
16        SUBROUTINE CG3D(        SUBROUTINE CG3D(
17       I                cg3d_b,       U                cg3d_b, cg3d_x,
18       U                cg3d_x,       O                firstResidual, lastResidual,
      O                firstResidual,  
      O                lastResidual,  
19       U                numIters,       U                numIters,
20       I                myIter, myThid )       I                myIter, myThid )
21  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
# Line 57  C     === Global data === Line 55  C     === Global data ===
55    
56  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
57  C     === Routine arguments ===  C     === Routine arguments ===
58  C     cg3d_b    :: The source term or "right hand side"  C     cg3d_b    :: The source term or "right hand side" (output: normalised RHS)
59  C     cg3d_x    :: The solution  C     cg3d_x    :: The solution (input: first guess)
60  C     firstResidual :: the initial residual before any iterations  C     firstResidual :: the initial residual before any iterations
61  C     lastResidual  :: the actual residual reached  C     minResidualSq :: the lowest residual reached (squared)
62  C     numIters  :: Entry: the maximum number of iterations allowed  CC    lastResidual  :: the actual residual reached
63  C               :: Exit:  the actual number of iterations used  C     numIters  :: Inp: the maximum number of iterations allowed
64    C                  Out: the actual number of iterations used
65    CC    nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
66    CC                 Out: iteration number corresponding to lowest residual
67  C     myIter    :: Current iteration number in simulation  C     myIter    :: Current iteration number in simulation
68  C     myThid    :: my Thread Id number  C     myThid    :: my Thread Id number
69        _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
# Line 77  C     myThid    :: my Thread Id number Line 78  C     myThid    :: my Thread Id number
78    
79  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
80  C     === Local variables ====  C     === Local variables ====
81  C     actualIts      :: Number of iterations taken  C     bi, bj     :: tile index in X and Y.
 C     actualResidual :: residual  
 C     bi,bj      :: tile indices  
 C     eta_qrN    :: Used in computing search directions  
 C     eta_qrNM1     suffix N and NM1 denote current and  
 C     cgBeta        previous iterations respectively.  
 C     alpha  
 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a  
 C                   useful debuggin/trouble shooting diagnostic.  
 C                   For neumann problems sumRHS needs to be ~0.  
 C                   or they converge at a non-zero residual.  
 C     err        :: Measure of residual of Ax - b, usually the norm.  
82  C     i, j, k    :: Loop counters  C     i, j, k    :: Loop counters
83  C     it3d       :: Loop counter for CG iterations  C     it3d       :: Loop counter for CG iterations
84    C     actualIts  :: actual CG iteration number
85    C     err_sq     :: Measure of the square of the residual of Ax - b.
86    C     eta_qrN    :: Used in computing search directions; suffix N and NM1
87    C     eta_qrNM1     denote current and previous iterations respectively.
88    C     cgBeta     :: coeff used to update conjugate direction vector "s".
89    C     alpha      :: coeff used to update solution & residual
90    C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
91    C                   debugging/trouble shooting diagnostic. For neumann problems
92    C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
93    C     cg2d_min   :: used to store solution corresponding to lowest residual.
94  C     msgBuf     :: Informational/error message buffer  C     msgBuf     :: Informational/error message buffer
       INTEGER actualIts  
       _RL     actualResidual  
95        INTEGER bi, bj        INTEGER bi, bj
96        INTEGER i, j, k, it3d        INTEGER i, j, k, it3d
97          INTEGER actualIts
98        INTEGER km1, kp1        INTEGER km1, kp1
99        _RL     maskM1, maskP1        _RL     maskM1, maskP1
100        _RL     err,    errTile(nSx,nSy)        _RL     cg3dTolerance_sq
101        _RL     eta_qrN,eta_qrNtile(nSx,nSy)        _RL     err_sq,  errTile(nSx,nSy)
102          _RL     eta_qrN, eta_qrNtile(nSx,nSy)
103        _RL     eta_qrNM1        _RL     eta_qrNM1
104        _RL     cgBeta        _RL     cgBeta
105        _RL     alpha , alphaTile(nSx,nSy)        _RL     alpha ,  alphaTile(nSx,nSy)
106        _RL     sumRHS, sumRHStile(nSx,nSy)        _RL     sumRHS,  sumRHStile(nSx,nSy)
107        _RL     rhsMax        _RL     rhsMax
108        _RL     rhsNorm        _RL     rhsNorm
109        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
# Line 115  C     msgBuf     :: Informational/error Line 115  C     msgBuf     :: Informational/error
115  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
116  CEOP  CEOP
117    
118    C--   Initialise auxiliary constant, some output variable
119          cg3dTolerance_sq = cg3dTargetResidual*cg3dTargetResidual
120        IF ( select_rStar .NE. 0 ) THEN        IF ( select_rStar .NE. 0 ) THEN
121          surfFac = freeSurfFac          surfFac = freeSurfFac
122        ELSE        ELSE
# Line 163  C--   Normalise RHS Line 165  C--   Normalise RHS
165        ENDDO        ENDDO
166    
167  C--   Update overlaps  C--   Update overlaps
 c     _EXCH_XYZ_RL( cg3d_b, myThid )  
168        _EXCH_XYZ_RL( cg3d_x, myThid )        _EXCH_XYZ_RL( cg3d_x, myThid )
169    
170  C--   Initial residual calculation (with free-Surface term)  C--   Initial residual calculation (with free-Surface term)
       err    = 0. _d 0  
       sumRHS = 0. _d 0  
171        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
172         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
173          errTile(bi,bj)    = 0. _d 0          errTile(bi,bj)    = 0. _d 0
# Line 229  C--   Initial residual calculation (with Line 228  C--   Initial residual calculation (with
228             sumRHStile(bi,bj) = sumRHStile(bi,bj)+cg3d_b(i,j,k,bi,bj)             sumRHStile(bi,bj) = sumRHStile(bi,bj)+cg3d_b(i,j,k,bi,bj)
229            ENDDO            ENDDO
230           ENDDO           ENDDO
231           DO J=1-1,sNy+1           DO j=0,sNy+1
232            DO I=1-1,sNx+1            DO i=0,sNx+1
233             cg3d_s(i,j,k,bi,bj) = 0.             cg3d_s(i,j,k,bi,bj) = 0.
234            ENDDO            ENDDO
235           ENDDO           ENDDO
# Line 238  C--   Initial residual calculation (with Line 237  C--   Initial residual calculation (with
237         ENDDO         ENDDO
238        ENDDO        ENDDO
239         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
 c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )  
240         CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )         CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
241         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
242        IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN        IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
243          CALL WRITE_FLD_S3D_RL(          CALL WRITE_FLD_S3D_RL(
244       I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )       I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )
245        ENDIF        ENDIF
246    
247          actualIts = 0
248          firstResidual = SQRT(err_sq)
249    
250        printResidual = .FALSE.        printResidual = .FALSE.
251        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
252          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
# Line 255  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myT Line 256  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myT
256          _END_MASTER( myThid )          _END_MASTER( myThid )
257        ENDIF        ENDIF
258    
259        actualIts      = 0        IF ( err_sq .LT. cg3dTolerance_sq ) GOTO 11
       actualResidual = SQRT(err)  
       firstResidual=actualResidual  
260    
261  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
262        DO 10 it3d=1, numIters        DO 10 it3d=1, numIters
263    
        IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11  
264  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
265  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
266  C      Note. On the next two loops over all tiles the inner loop ranges  C      Note. On the next two loops over all tiles the inner loop ranges
267  C            in sNx and sNy are expanded by 1 to avoid a communication  C            in sNx and sNy are expanded by 1 to avoid a communication
268  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
269  C            want eta_qrN for the interior points.  C            want eta_qrN for the interior points.
        eta_qrN = 0. _d 0  
270         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
271          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
272           eta_qrNtile(bi,bj) = 0. _d 0           eta_qrNtile(bi,bj) = 0. _d 0
# Line 277  C            want eta_qrN for the interi Line 274  C            want eta_qrN for the interi
274  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
275  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
276  #endif /* TARGET_NEC_SX */  #endif /* TARGET_NEC_SX */
277            DO j=1-1,sNy+1            DO j=0,sNy+1
278             DO i=1-1,sNx+1             DO i=0,sNx+1
279              cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)              cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
280       &                        *cg3d_r(i,j,k,bi,bj)       &                        *cg3d_r(i,j,k,bi,bj)
281             ENDDO             ENDDO
# Line 288  C            want eta_qrN for the interi Line 285  C            want eta_qrN for the interi
285  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
286  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
287  #endif /* TARGET_NEC_SX */  #endif /* TARGET_NEC_SX */
288            DO j=1-1,sNy+1            DO j=0,sNy+1
289             DO i=1-1,sNx+1             DO i=0,sNx+1
290              cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)              cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
291       &                      *( cg3d_r(i,j,k,bi,bj)       &                      *( cg3d_r(i,j,k,bi,bj)
292       &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)       &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)
# Line 312  C            want eta_qrN for the interi Line 309  C            want eta_qrN for the interi
309  #ifdef TARGET_NEC_SX  #ifdef TARGET_NEC_SX
310  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
311  #endif /* TARGET_NEC_SX */  #endif /* TARGET_NEC_SX */
312            DO j=1-1,sNy+1            DO j=0,sNy+1
313             DO i=1-1,sNx+1             DO i=0,sNx+1
314              cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)              cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
315       &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)       &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)
316             ENDDO             ENDDO
# Line 332  C            want eta_qrN for the interi Line 329  C            want eta_qrN for the interi
329         ENDDO         ENDDO
330    
331         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  
 CcnhDebugEnds  
332         cgBeta   = eta_qrN/eta_qrNM1         cgBeta   = eta_qrN/eta_qrNM1
333  CcnhDebugStarts  CcnhDebugStarts
334  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta  c      WRITE(*,*) ' CG3D: Iteration ', it3d-1,
335    c    &            ' eta_qrN=', eta_qrN, ' beta=', cgBeta
336  CcnhDebugEnds  CcnhDebugEnds
337         eta_qrNM1 = eta_qrN         eta_qrNM1 = eta_qrN
338    
339         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
340          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
341           DO k=1,Nr           DO k=1,Nr
342            DO j=1-1,sNy+1            DO j=0,sNy+1
343             DO i=1-1,sNx+1             DO i=0,sNx+1
344              cg3d_s(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)              cg3d_s(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
345       &                   + cgBeta*cg3d_s(i,j,k,bi,bj)       &                   + cgBeta*cg3d_s(i,j,k,bi,bj)
346             ENDDO             ENDDO
# Line 356  CcnhDebugEnds Line 351  CcnhDebugEnds
351    
352  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
353  C==    q = A.s  C==    q = A.s
        alpha = 0. _d 0  
354         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
355          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
356           alphaTile(bi,bj) = 0. _d 0           alphaTile(bi,bj) = 0. _d 0
# Line 476  C==    q = A.s Line 470  C==    q = A.s
470         ENDDO         ENDDO
471         CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )         CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
472  CcnhDebugStarts  CcnhDebugStarts
473  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  c      WRITE(*,*) ' CG3D: Iteration ', it3d-1,
474    c    &            ' SUM(s*q)=', alpha, ' alpha=', eta_qrN/alpha
475  CcnhDebugEnds  CcnhDebugEnds
476         alpha = eta_qrN/alpha         alpha = eta_qrN/alpha
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  
 CcnhDebugEnds  
477    
478  C==    Update solution and residual vectors  C==    Update simultaneously solution and residual vectors (and Iter number)
479  C      Now compute "interior" points.  C      Now compute "interior" points.
        err = 0. _d 0  
480         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
481          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
482           errTile(bi,bj) = 0. _d 0           errTile(bi,bj) = 0. _d 0
# Line 506  C      Now compute "interior" points. Line 497  C      Now compute "interior" points.
497           ENDDO           ENDDO
498          ENDDO          ENDDO
499         ENDDO         ENDDO
500           actualIts = it3d
501    
502         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
        err = SQRT(err)  
        actualIts      = it3d  
        actualResidual = err  
503         IF ( printResidual ) THEN         IF ( printResidual ) THEN
504          IF ( MOD( it3d-1, printResidualFreq ).EQ.0 ) THEN          IF ( MOD( it3d-1, printResidualFreq ).EQ.0 ) THEN
505           WRITE(msgBuf,'(A,I6,A,1PE21.14)')           WRITE(msgBuf,'(A,I6,A,1PE21.14)')
506       &    ' cg3d: iter=', actualIts, ' ; resid.= ', actualResidual       &    ' cg3d: iter=', it3d, ' ; resid.= ', SQRT(err_sq)
507           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
508       &                       SQUEEZE_RIGHT, myThid )       &                       SQUEEZE_RIGHT, myThid )
509          ENDIF          ENDIF
510         ENDIF         ENDIF
511         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( err_sq .LT. cg3dTolerance_sq ) GOTO 11
512         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
513    
514     10 CONTINUE     10 CONTINUE
# Line 543  C--   Un-normalise the answer Line 532  C--   Un-normalise the answer
532         ENDDO         ENDDO
533        ENDDO        ENDDO
534    
535        lastResidual = actualResidual  C--   Return parameters to caller
536          lastResidual = SQRT(err_sq)
537        numIters = actualIts        numIters = actualIts
538    
539  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22