Back to home page

MITgcm
 
 

    

File indexing completed on 2021-11-09 06:14:42 UTC

view on githubraw file Latest commit d68b36a4 on 2021-11-08 15:58:39 UTC
23bce0bbb8 Mart*0001 #include "CPP_OPTIONS.h"
56f3d2cc76 Mart*0002 #ifdef TARGET_NEC_SX
                0003 C     set a sensible default for the outer loop unrolling parameter that can
                0004 C     be overriden in the Makefile with the DEFINES macro or in CPP_OPTIONS.h
                0005 #ifndef CG2D_OUTERLOOPITERS
                0006 # define CG2D_OUTERLOOPITERS 10
                0007 #endif
                0008 #endif /* TARGET_NEC_SX */
                0009 
23bce0bbb8 Mart*0010 CBOP
                0011 C     !ROUTINE: CG2D_SR
                0012 C     !INTERFACE:
                0013       SUBROUTINE CG2D_SR(
c90853988b Jean*0014      U                cg2d_b, cg2d_x,
                0015      O                firstResidual, minResidualSq, lastResidual,
                0016      U                numIters, nIterMin,
23bce0bbb8 Mart*0017      I                myThid )
                0018 C     !DESCRIPTION: \bv
aecc8b0f47 Mart*0019 C     *================================================================*
                0020 C     | SUBROUTINE CG2D_SR
                0021 C     | o Two-dimensional grid problem conjugate-gradient inverter
                0022 C     |   (with preconditioner).
                0023 C     | o This version used the more efficient (with less global
                0024 C     |   communications) "Single Reduction" (SR) implementation from
                0025 C     |   d Azevedo, Eijkhout, and Romine (Lapack Working Note 56, 1999)
                0026 C     | Author:  C. Wolfe, November 2009, clwolfe@ucsd.edu
                0027 C     *================================================================*
23bce0bbb8 Mart*0028 C     | Con. grad is an iterative procedure for solving Ax = b.
                0029 C     | It requires the A be symmetric.
aecc8b0f47 Mart*0030 C     | This implementation assumes A is a five-diagonal matrix
                0031 C     | of the form that arises in the discrete representation of
                0032 C     | the del^2 operator in a two-dimensional space.
                0033 C     *================================================================*
23bce0bbb8 Mart*0034 C     \ev
                0035 
                0036 C     !USES:
                0037       IMPLICIT NONE
                0038 C     === Global data ===
                0039 #include "SIZE.h"
                0040 #include "EEPARAMS.h"
                0041 #include "PARAMS.h"
                0042 #include "CG2D.h"
                0043 
                0044 C     !INPUT/OUTPUT PARAMETERS:
c90853988b Jean*0045 C     cg2d_b    :: The source term or "right hand side" (output: normalised RHS)
                0046 C     cg2d_x    :: The solution (input: first guess)
23bce0bbb8 Mart*0047 C     firstResidual :: the initial residual before any iterations
c90853988b Jean*0048 C     minResidualSq :: the lowest residual reached (squared)
23bce0bbb8 Mart*0049 C     lastResidual  :: the actual residual reached
c90853988b Jean*0050 C     numIters  :: Inp: the maximum number of iterations allowed
                0051 C                  Out: the actual number of iterations used
                0052 C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
                0053 C                  Out: iteration number corresponding to lowest residual
                0054 C     myThid    :: Thread on which I am working.
23bce0bbb8 Mart*0055       _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0056       _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0057       _RL  firstResidual
c90853988b Jean*0058       _RL  minResidualSq
23bce0bbb8 Mart*0059       _RL  lastResidual
                0060       INTEGER numIters
c90853988b Jean*0061       INTEGER nIterMin
23bce0bbb8 Mart*0062       INTEGER myThid
                0063 
                0064 #ifdef ALLOW_SRCG
aecc8b0f47 Mart*0065 C--   Local variables in COMMON block (shared by all threads)
                0066 C     sumPhi     :: needed to call global_vec_sum_r8 when mutli-threaded
                0067       _RL    sumPhi(3,nSx,nSy)
                0068       COMMON /CG2D_SR_LOCAL/ sumPhi
23bce0bbb8 Mart*0069 
                0070 C     !LOCAL VARIABLES:
c90853988b Jean*0071 C     bi, bj     :: tile index in X and Y.
                0072 C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
                0073 C     actualIts  :: actual CG iteration number
                0074 C     err_sq     :: Measure of the square of the residual of Ax - b.
                0075 C     eta_qrN    :: Used in computing search directions; suffix N and NM1
                0076 C     eta_qrNM1     denote current and previous iterations respectively.
                0077 C     cgBeta     :: coeff used to update conjugate direction vector "s".
                0078 C     alpha      :: coeff used to update solution & residual
                0079 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
                0080 C                   debugging/trouble shooting diagnostic. For neumann problems
                0081 C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
                0082 C     cg2d_min   :: used to store solution corresponding to lowest residual.
                0083 C     msgBuf     :: Informational/error message buffer
aecc8b0f47 Mart*0084 C--   local working array (used to be in CG2D.h common block:
                0085 C     cg2d_y     :: residual scaled by preconditioner
                0086 C     cg2d_v     :: z times the operator
                0087 C     cg2d_q     :: Intermediate matrix-vector product term
                0088 C     cg2d_r     ::   *same*
                0089 C     cg2d_s     ::   *same*
23bce0bbb8 Mart*0090       INTEGER bi, bj
c90853988b Jean*0091       INTEGER i, j, it2d
                0092 c     INTEGER actualIts
                0093       _RL    cg2dTolerance_sq
                0094       _RL    err_sq,  errTile(nSx,nSy)
                0095       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
23bce0bbb8 Mart*0096       _RL    eta_qrNM1
                0097       _RL    cgBeta
                0098       _RL    alpha,  alphaTile(nSx,nSy)
e1fb02e8f0 Jean*0099       _RL    sigma
23bce0bbb8 Mart*0100       _RL    delta,  deltaTile(nSx,nSy)
                0101       _RL    sumRHS, sumRHStile(nSx,nSy)
                0102       _RL    rhsMax
                0103       _RL    rhsNorm
c90853988b Jean*0104       _RL    cg2d_min(1:sNx,1:sNy,nSx,nSy)
d68b36a485 Mart*0105       _RL    cg2d_v(1:sNx,1:sNy,nSx,nSy)
                0106       _RL    cg2d_q(1:sNx,1:sNy,nSx,nSy)
aecc8b0f47 Mart*0107       _RL    cg2d_y(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0108       _RL    cg2d_r(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0109       _RL    cg2d_s(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
23bce0bbb8 Mart*0110       CHARACTER*(MAX_LEN_MBUF) msgBuf
764aea3440 Jean*0111       LOGICAL printResidual
23bce0bbb8 Mart*0112 CEOP
                0113 
c90853988b Jean*0114 C--   Initialise auxiliary constant, some output variable and inverter
                0115       cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
                0116       minResidualSq = -1. _d 0
                0117       eta_qrNM1     =  1. _d 0
23bce0bbb8 Mart*0118 
                0119 C--   Normalise RHS
                0120       rhsMax = 0. _d 0
                0121       DO bj=myByLo(myThid),myByHi(myThid)
                0122        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0123         DO j=1,sNy
                0124          DO i=1,sNx
                0125           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0126           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
23bce0bbb8 Mart*0127          ENDDO
                0128         ENDDO
                0129        ENDDO
                0130       ENDDO
                0131 
                0132       IF (cg2dNormaliseRHS) THEN
                0133 C-  Normalise RHS :
                0134       _GLOBAL_MAX_RL( rhsMax, myThid )
                0135       rhsNorm = 1. _d 0
                0136       IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0137       DO bj=myByLo(myThid),myByHi(myThid)
                0138        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0139         DO j=1,sNy
                0140          DO i=1,sNx
                0141           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0142           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
23bce0bbb8 Mart*0143          ENDDO
                0144         ENDDO
                0145        ENDDO
                0146       ENDDO
                0147 C- end Normalise RHS
                0148       ENDIF
                0149 
                0150 C--   Update overlaps
                0151       CALL EXCH_XY_RL( cg2d_x, myThid )
                0152 
                0153 C--   Initial residual calculation
                0154       DO bj=myByLo(myThid),myByHi(myThid)
                0155        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0156 C-    Initialise local working arrays:
d68b36a485 Mart*0157         DO j=1,sNy
                0158          DO i=1,sNx
                0159           cg2d_v(i,j,bi,bj) = 0. _d 0
                0160           cg2d_q(i,j,bi,bj) = 0. _d 0
                0161          ENDDO
                0162         ENDDO
ff51fd05c6 Jean*0163         DO j=0,sNy+1
                0164          DO i=0,sNx+1
                0165           cg2d_y(i,j,bi,bj) = 0. _d 0
                0166           cg2d_r(i,j,bi,bj) = 0. _d 0
                0167           cg2d_s(i,j,bi,bj) = 0. _d 0
                0168          ENDDO
                0169         ENDDO
c90853988b Jean*0170         IF ( nIterMin.GE.0 ) THEN
ff51fd05c6 Jean*0171 C-    Initialise saved solution
c90853988b Jean*0172          DO j=1,sNy
                0173           DO i=1,sNx
ff51fd05c6 Jean*0174            cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
c90853988b Jean*0175           ENDDO
                0176          ENDDO
                0177         ENDIF
23bce0bbb8 Mart*0178         sumRHStile(bi,bj) = 0. _d 0
                0179         errTile(bi,bj)    = 0. _d 0
56f3d2cc76 Mart*0180 #ifdef TARGET_NEC_SX
                0181 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0182 #endif /* TARGET_NEC_SX */
c90853988b Jean*0183         DO j=1,sNy
                0184          DO i=1,sNx
                0185           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0186      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0187      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0188      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0189      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0190      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0191      &    )
                0192           errTile(bi,bj)    = errTile(bi,bj)
c90853988b Jean*0193      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0194           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
23bce0bbb8 Mart*0195          ENDDO
                0196         ENDDO
                0197        ENDDO
                0198       ENDDO
ff02675122 Jean*0199 
23bce0bbb8 Mart*0200       CALL EXCH_S3D_RL( cg2d_r, 1,  myThid )
ff02675122 Jean*0201 
23bce0bbb8 Mart*0202       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
c90853988b Jean*0203       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
                0204 
                0205       it2d      = 0
                0206 c     actualIts = 0
                0207       firstResidual = SQRT(err_sq)
                0208       IF ( nIterMin.GE.0 ) THEN
                0209         nIterMin = 0
                0210         minResidualSq = err_sq
                0211       ENDIF
23bce0bbb8 Mart*0212 
764aea3440 Jean*0213       printResidual = .FALSE.
23bce0bbb8 Mart*0214       IF ( debugLevel .GE. debLevZero ) THEN
                0215         _BEGIN_MASTER( myThid )
764aea3440 Jean*0216         printResidual = printResidualFreq.GE.1
23bce0bbb8 Mart*0217         WRITE(standardmessageunit,'(A,1P2E22.14)')
                0218      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
                0219         _END_MASTER( myThid )
                0220       ENDIF
ff02675122 Jean*0221 
c90853988b Jean*0222       IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
ff02675122 Jean*0223 
23bce0bbb8 Mart*0224 C DER (1999) do one iteration outside of the loop to start things up.
                0225 C--    Solve preconditioning equation and update
                0226 C--    conjugate direction vector "s".
                0227       DO bj=myByLo(myThid),myByHi(myThid)
                0228        DO bi=myBxLo(myThid),myBxHi(myThid)
                0229         eta_qrNtile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0230 #ifdef TARGET_NEC_SX
                0231 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0232 #endif /* TARGET_NEC_SX */
c90853988b Jean*0233         DO j=1,sNy
                0234          DO i=1,sNx
                0235            cg2d_y(i,j,bi,bj) =
                0236      &       pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0237      &      +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0238      &      +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0239      &      +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0240      &      +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
                0241            cg2d_s(i,j,bi,bj)  = cg2d_y(i,j,bi,bj)
23bce0bbb8 Mart*0242            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
c90853988b Jean*0243      &      +cg2d_y(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0244          ENDDO
                0245         ENDDO
                0246        ENDDO
                0247       ENDDO
ff02675122 Jean*0248 
23bce0bbb8 Mart*0249       CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
                0250       CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
ff02675122 Jean*0251 
23bce0bbb8 Mart*0252       eta_qrNM1 = eta_qrN
                0253 
                0254 C==    Evaluate laplace operator on conjugate gradient vector
                0255 C==    q = A.s
                0256       DO bj=myByLo(myThid),myByHi(myThid)
                0257        DO bi=myBxLo(myThid),myBxHi(myThid)
                0258         alphaTile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0259 #ifdef TARGET_NEC_SX
                0260 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0261 #endif /* TARGET_NEC_SX */
c90853988b Jean*0262         DO j=1,sNy
                0263          DO i=1,sNx
                0264           cg2d_q(i,j,bi,bj) =
                0265      &      aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0266      &     +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0267      &     +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0268      &     +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0269      &     +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0270           alphaTile(bi,bj) = alphaTile(bi,bj)
c90853988b Jean*0271      &                     + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0272          ENDDO
                0273         ENDDO
                0274        ENDDO
                0275       ENDDO
ff02675122 Jean*0276 
23bce0bbb8 Mart*0277       CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
ff02675122 Jean*0278 
23bce0bbb8 Mart*0279       sigma = eta_qrN/alpha
                0280 
c90853988b Jean*0281 C==   Update simultaneously solution and residual vectors (and Iter number)
23bce0bbb8 Mart*0282 C     Now compute "interior" points.
                0283       DO bj=myByLo(myThid),myByHi(myThid)
                0284        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0285         DO j=1,sNy
                0286          DO i=1,sNx
                0287           cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+sigma*cg2d_s(i,j,bi,bj)
                0288           cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-sigma*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0289          ENDDO
                0290         ENDDO
                0291        ENDDO
                0292       ENDDO
c90853988b Jean*0293 c     actualIts = 1
23bce0bbb8 Mart*0294 
                0295       CALL EXCH_S3D_RL( cg2d_r,1, myThid )
                0296 
                0297 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c90853988b Jean*0298       DO 10 it2d=1, numIters-1
                0299 c      actualIts = it2d
23bce0bbb8 Mart*0300 
                0301 C--    Solve preconditioning equation and update
                0302 C--    conjugate direction vector "s".
                0303 C--    z = M^-1 r
                0304        DO bj=myByLo(myThid),myByHi(myThid)
                0305         DO bi=myBxLo(myThid),myBxHi(myThid)
56f3d2cc76 Mart*0306 #ifdef TARGET_NEC_SX
                0307 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0308 #endif /* TARGET_NEC_SX */
c90853988b Jean*0309          DO j=1,sNy
                0310           DO i=1,sNx
                0311            cg2d_y(i,j,bi,bj) =
                0312      &       pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0313      &      +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0314      &      +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0315      &      +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0316      &      +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
23bce0bbb8 Mart*0317           ENDDO
                0318          ENDDO
                0319         ENDDO
                0320        ENDDO
                0321 
                0322        CALL EXCH_S3D_RL( cg2d_y, 1, myThid )
ff02675122 Jean*0323 
23bce0bbb8 Mart*0324 C==    v = A.z
                0325 C--    eta_qr = <z,r>
                0326 C--    delta = <z,v>
                0327 C--    Do the error calcuation here to consolidate global reductions
                0328        DO bj=myByLo(myThid),myByHi(myThid)
                0329         DO bi=myBxLo(myThid),myBxHi(myThid)
                0330          eta_qrNtile(bi,bj) = 0. _d 0
                0331          deltaTile(bi,bj)   = 0. _d 0
                0332          errTile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0333 #ifdef TARGET_NEC_SX
                0334 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0335 #endif /* TARGET_NEC_SX */
c90853988b Jean*0336          DO j=1,sNy
                0337           DO i=1,sNx
                0338            cg2d_v(i,j,bi,bj) =
                0339      &       aW2d(i  ,j  ,bi,bj)*cg2d_y(i-1,j  ,bi,bj)
                0340      &      +aW2d(i+1,j  ,bi,bj)*cg2d_y(i+1,j  ,bi,bj)
                0341      &      +aS2d(i  ,j  ,bi,bj)*cg2d_y(i  ,j-1,bi,bj)
                0342      &      +aS2d(i  ,j+1,bi,bj)*cg2d_y(i  ,j+1,bi,bj)
                0343      &      +aC2d(i  ,j  ,bi,bj)*cg2d_y(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0344            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
c90853988b Jean*0345      &      +cg2d_y(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0346            deltaTile(bi,bj) = deltaTile(bi,bj)
c90853988b Jean*0347      &      +cg2d_y(i,j,bi,bj)*cg2d_v(i,j,bi,bj)
23bce0bbb8 Mart*0348            errTile(bi,bj) = errTile(bi,bj)
c90853988b Jean*0349      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0350           ENDDO
                0351          ENDDO
                0352          sumPhi(1,bi,bj) = eta_qrNtile(bi,bj)
                0353          sumPhi(2,bi,bj) = deltaTile(bi,bj)
                0354          sumPhi(3,bi,bj) = errTile(bi,bj)
                0355         ENDDO
ff02675122 Jean*0356        ENDDO
                0357 
c90853988b Jean*0358 C      CALL GLOBAL_SUM_TILE_RL( eta_qrNtile, eta_qrN,myThid )
                0359 C      CALL GLOBAL_SUM_TILE_RL( deltaTile, delta, myThid )
                0360 C      CALL GLOBAL_SUM_TILE_RL( errTile,  err_sq, myThid )
                0361        CALL GLOBAL_VEC_SUM_R8( 3, 3, sumPhi, myThid )
ff02675122 Jean*0362 
23bce0bbb8 Mart*0363        eta_qrN = sumPhi(1,1,1)
                0364        delta   = sumPhi(2,1,1)
c90853988b Jean*0365        err_sq  = sumPhi(3,1,1)
ff02675122 Jean*0366 
764aea3440 Jean*0367        IF ( printResidual ) THEN
                0368         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0369          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
c90853988b Jean*0370      &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
764aea3440 Jean*0371          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0372      &                       SQUEEZE_RIGHT, myThid )
                0373         ENDIF
23bce0bbb8 Mart*0374        ENDIF
c90853988b Jean*0375        IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
                0376        IF ( err_sq .LT. minResidualSq ) THEN
                0377 C-     Store lowest residual solution
                0378          minResidualSq = err_sq
                0379          nIterMin = it2d
                0380          DO bj=myByLo(myThid),myByHi(myThid)
                0381           DO bi=myBxLo(myThid),myBxHi(myThid)
                0382            DO j=1,sNy
                0383             DO i=1,sNx
                0384              cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0385             ENDDO
                0386            ENDDO
                0387           ENDDO
                0388          ENDDO
                0389        ENDIF
ff02675122 Jean*0390 
23bce0bbb8 Mart*0391        cgBeta   = eta_qrN/eta_qrNM1
                0392        eta_qrNM1 = eta_qrN
c90853988b Jean*0393        alpha = delta - (cgBeta*cgBeta)*alpha
23bce0bbb8 Mart*0394        sigma = eta_qrN/alpha
ff02675122 Jean*0395 
c90853988b Jean*0396 C==    Update simultaneously solution and residual vectors (and Iter number)
23bce0bbb8 Mart*0397        DO bj=myByLo(myThid),myByHi(myThid)
                0398         DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0399          DO j=1,sNy
                0400           DO i=1,sNx
                0401            cg2d_s(i,j,bi,bj) = cg2d_y(i,j,bi,bj)
                0402      &                       + cgBeta*cg2d_s(i,j,bi,bj)
                0403            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0404      &                       + sigma*cg2d_s(i,j,bi,bj)
                0405            cg2d_q(i,j,bi,bj) = cg2d_v(i,j,bi,bj)
                0406      &                       + cgBeta*cg2d_q(i,j,bi,bj)
                0407            cg2d_r(i,j,bi,bj) = cg2d_r(i,j,bi,bj)
                0408      &                       - sigma*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0409           ENDDO
                0410          ENDDO
                0411         ENDDO
                0412        ENDDO
c90853988b Jean*0413 c      actualIts = it2d + 1
23bce0bbb8 Mart*0414 
                0415        CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0416 
                0417    10 CONTINUE
c90853988b Jean*0418       DO bj=myByLo(myThid),myByHi(myThid)
                0419        DO bi=myBxLo(myThid),myBxHi(myThid)
                0420         errTile(bi,bj) = 0. _d 0
                0421         DO j=1,sNy
                0422          DO i=1,sNx
                0423           errTile(bi,bj) = errTile(bi,bj)
                0424      &                   + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0425          ENDDO
                0426         ENDDO
                0427        ENDDO
                0428       ENDDO
                0429       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
23bce0bbb8 Mart*0430    11 CONTINUE
                0431 
c90853988b Jean*0432       IF ( nIterMin.GE.0 .AND. err_sq .GT. minResidualSq ) THEN
                0433 C-    use the lowest residual solution (instead of current one = last residual)
                0434         DO bj=myByLo(myThid),myByHi(myThid)
                0435          DO bi=myBxLo(myThid),myBxHi(myThid)
                0436           DO j=1,sNy
                0437            DO i=1,sNx
                0438              cg2d_x(i,j,bi,bj) = cg2d_min(i,j,bi,bj)
                0439            ENDDO
                0440           ENDDO
                0441          ENDDO
                0442         ENDDO
                0443       ENDIF
                0444 
23bce0bbb8 Mart*0445       IF (cg2dNormaliseRHS) THEN
                0446 C--   Un-normalise the answer
                0447         DO bj=myByLo(myThid),myByHi(myThid)
                0448          DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0449           DO j=1,sNy
                0450            DO i=1,sNx
                0451             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
23bce0bbb8 Mart*0452            ENDDO
                0453           ENDDO
                0454          ENDDO
                0455         ENDDO
                0456       ENDIF
                0457 
                0458 C--   Return parameters to caller
c90853988b Jean*0459       lastResidual = SQRT(err_sq)
                0460       numIters = it2d
                0461 c     numIters = actualIts
23bce0bbb8 Mart*0462 
                0463 #endif /* ALLOW_SRCG */
                0464       RETURN
                0465       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!