Back to home page

MITgcm
 
 

    

File indexing completed on 2024-09-20 05:10:17 UTC

view on githubraw file Latest commit e6e223b2 on 2024-09-19 22:01:15 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
                0065 C     !LOCAL VARIABLES:
c90853988b Jean*0066 C     bi, bj     :: tile index in X and Y.
                0067 C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
                0068 C     actualIts  :: actual CG iteration number
                0069 C     err_sq     :: Measure of the square of the residual of Ax - b.
                0070 C     eta_qrN    :: Used in computing search directions; suffix N and NM1
                0071 C     eta_qrNM1     denote current and previous iterations respectively.
                0072 C     cgBeta     :: coeff used to update conjugate direction vector "s".
                0073 C     alpha      :: coeff used to update solution & residual
                0074 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
                0075 C                   debugging/trouble shooting diagnostic. For neumann problems
                0076 C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
5237154b93 Jean*0077 C     sumVec     :: output vector for GLOBAL_SUM_VECTOR call
                0078 C     tiledVec   :: input  vector for GLOBAL_SUM_VECTOR call
c90853988b Jean*0079 C     cg2d_min   :: used to store solution corresponding to lowest residual.
                0080 C     msgBuf     :: Informational/error message buffer
aecc8b0f47 Mart*0081 C--   local working array (used to be in CG2D.h common block:
                0082 C     cg2d_y     :: residual scaled by preconditioner
                0083 C     cg2d_v     :: z times the operator
                0084 C     cg2d_q     :: Intermediate matrix-vector product term
                0085 C     cg2d_r     ::   *same*
                0086 C     cg2d_s     ::   *same*
23bce0bbb8 Mart*0087       INTEGER bi, bj
c90853988b Jean*0088       INTEGER i, j, it2d
                0089 c     INTEGER actualIts
                0090       _RL    err_sq,  errTile(nSx,nSy)
                0091       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
23bce0bbb8 Mart*0092       _RL    eta_qrNM1
                0093       _RL    cgBeta
                0094       _RL    alpha,  alphaTile(nSx,nSy)
e1fb02e8f0 Jean*0095       _RL    sigma
23bce0bbb8 Mart*0096       _RL    delta,  deltaTile(nSx,nSy)
                0097       _RL    sumRHS, sumRHStile(nSx,nSy)
5237154b93 Jean*0098       _RL    sumVec(3), tiledVec(nSx,nSy,3)
23bce0bbb8 Mart*0099       _RL    rhsMax
                0100       _RL    rhsNorm
c90853988b Jean*0101       _RL    cg2d_min(1:sNx,1:sNy,nSx,nSy)
d68b36a485 Mart*0102       _RL    cg2d_v(1:sNx,1:sNy,nSx,nSy)
                0103       _RL    cg2d_q(1:sNx,1:sNy,nSx,nSy)
aecc8b0f47 Mart*0104       _RL    cg2d_y(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0105       _RL    cg2d_r(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0106       _RL    cg2d_s(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
23bce0bbb8 Mart*0107       CHARACTER*(MAX_LEN_MBUF) msgBuf
764aea3440 Jean*0108       LOGICAL printResidual
23bce0bbb8 Mart*0109 CEOP
                0110 
c90853988b Jean*0111 C--   Initialise auxiliary constant, some output variable and inverter
                0112       minResidualSq = -1. _d 0
                0113       eta_qrNM1     =  1. _d 0
23bce0bbb8 Mart*0114 
                0115 C--   Normalise RHS
                0116       rhsMax = 0. _d 0
                0117       DO bj=myByLo(myThid),myByHi(myThid)
                0118        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0119         DO j=1,sNy
                0120          DO i=1,sNx
                0121           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0122           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
23bce0bbb8 Mart*0123          ENDDO
                0124         ENDDO
                0125        ENDDO
                0126       ENDDO
                0127 
                0128       IF (cg2dNormaliseRHS) THEN
                0129 C-  Normalise RHS :
                0130       _GLOBAL_MAX_RL( rhsMax, myThid )
                0131       rhsNorm = 1. _d 0
                0132       IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0133       DO bj=myByLo(myThid),myByHi(myThid)
                0134        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0135         DO j=1,sNy
                0136          DO i=1,sNx
                0137           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0138           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
23bce0bbb8 Mart*0139          ENDDO
                0140         ENDDO
                0141        ENDDO
                0142       ENDDO
                0143 C- end Normalise RHS
                0144       ENDIF
                0145 
                0146 C--   Update overlaps
                0147       CALL EXCH_XY_RL( cg2d_x, myThid )
                0148 
                0149 C--   Initial residual calculation
                0150       DO bj=myByLo(myThid),myByHi(myThid)
                0151        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0152 C-    Initialise local working arrays:
d68b36a485 Mart*0153         DO j=1,sNy
                0154          DO i=1,sNx
                0155           cg2d_v(i,j,bi,bj) = 0. _d 0
                0156           cg2d_q(i,j,bi,bj) = 0. _d 0
                0157          ENDDO
                0158         ENDDO
ff51fd05c6 Jean*0159         DO j=0,sNy+1
                0160          DO i=0,sNx+1
                0161           cg2d_y(i,j,bi,bj) = 0. _d 0
                0162           cg2d_r(i,j,bi,bj) = 0. _d 0
                0163           cg2d_s(i,j,bi,bj) = 0. _d 0
                0164          ENDDO
                0165         ENDDO
c90853988b Jean*0166         IF ( nIterMin.GE.0 ) THEN
ff51fd05c6 Jean*0167 C-    Initialise saved solution
c90853988b Jean*0168          DO j=1,sNy
                0169           DO i=1,sNx
ff51fd05c6 Jean*0170            cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
c90853988b Jean*0171           ENDDO
                0172          ENDDO
                0173         ENDIF
23bce0bbb8 Mart*0174         sumRHStile(bi,bj) = 0. _d 0
                0175         errTile(bi,bj)    = 0. _d 0
56f3d2cc76 Mart*0176 #ifdef TARGET_NEC_SX
                0177 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0178 #endif /* TARGET_NEC_SX */
c90853988b Jean*0179         DO j=1,sNy
                0180          DO i=1,sNx
                0181           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0182      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0183      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0184      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0185      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0186      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0187      &    )
                0188           errTile(bi,bj)    = errTile(bi,bj)
c90853988b Jean*0189      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0190           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
23bce0bbb8 Mart*0191          ENDDO
                0192         ENDDO
                0193        ENDDO
                0194       ENDDO
ff02675122 Jean*0195 
23bce0bbb8 Mart*0196       CALL EXCH_S3D_RL( cg2d_r, 1,  myThid )
ff02675122 Jean*0197 
23bce0bbb8 Mart*0198       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
c90853988b Jean*0199       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
                0200 
                0201       it2d      = 0
                0202 c     actualIts = 0
                0203       firstResidual = SQRT(err_sq)
                0204       IF ( nIterMin.GE.0 ) THEN
                0205         nIterMin = 0
                0206         minResidualSq = err_sq
                0207       ENDIF
23bce0bbb8 Mart*0208 
764aea3440 Jean*0209       printResidual = .FALSE.
23bce0bbb8 Mart*0210       IF ( debugLevel .GE. debLevZero ) THEN
                0211         _BEGIN_MASTER( myThid )
764aea3440 Jean*0212         printResidual = printResidualFreq.GE.1
5237154b93 Jean*0213         WRITE(standardMessageUnit,'(A,1P2E22.14)')
23bce0bbb8 Mart*0214      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
                0215         _END_MASTER( myThid )
                0216       ENDIF
ff02675122 Jean*0217 
c90853988b Jean*0218       IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
ff02675122 Jean*0219 
23bce0bbb8 Mart*0220 C DER (1999) do one iteration outside of the loop to start things up.
                0221 C--    Solve preconditioning equation and update
                0222 C--    conjugate direction vector "s".
                0223       DO bj=myByLo(myThid),myByHi(myThid)
                0224        DO bi=myBxLo(myThid),myBxHi(myThid)
                0225         eta_qrNtile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0226 #ifdef TARGET_NEC_SX
                0227 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0228 #endif /* TARGET_NEC_SX */
c90853988b Jean*0229         DO j=1,sNy
                0230          DO i=1,sNx
                0231            cg2d_y(i,j,bi,bj) =
                0232      &       pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0233      &      +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0234      &      +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0235      &      +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0236      &      +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
                0237            cg2d_s(i,j,bi,bj)  = cg2d_y(i,j,bi,bj)
23bce0bbb8 Mart*0238            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
c90853988b Jean*0239      &      +cg2d_y(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0240          ENDDO
                0241         ENDDO
                0242        ENDDO
                0243       ENDDO
ff02675122 Jean*0244 
23bce0bbb8 Mart*0245       CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
                0246       CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
ff02675122 Jean*0247 
23bce0bbb8 Mart*0248       eta_qrNM1 = eta_qrN
                0249 
                0250 C==    Evaluate laplace operator on conjugate gradient vector
                0251 C==    q = A.s
                0252       DO bj=myByLo(myThid),myByHi(myThid)
                0253        DO bi=myBxLo(myThid),myBxHi(myThid)
                0254         alphaTile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0255 #ifdef TARGET_NEC_SX
                0256 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0257 #endif /* TARGET_NEC_SX */
c90853988b Jean*0258         DO j=1,sNy
                0259          DO i=1,sNx
                0260           cg2d_q(i,j,bi,bj) =
                0261      &      aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0262      &     +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0263      &     +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0264      &     +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0265      &     +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0266           alphaTile(bi,bj) = alphaTile(bi,bj)
c90853988b Jean*0267      &                     + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0268          ENDDO
                0269         ENDDO
                0270        ENDDO
                0271       ENDDO
ff02675122 Jean*0272 
23bce0bbb8 Mart*0273       CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
ff02675122 Jean*0274 
23bce0bbb8 Mart*0275       sigma = eta_qrN/alpha
                0276 
c90853988b Jean*0277 C==   Update simultaneously solution and residual vectors (and Iter number)
23bce0bbb8 Mart*0278 C     Now compute "interior" points.
                0279       DO bj=myByLo(myThid),myByHi(myThid)
                0280        DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0281         DO j=1,sNy
                0282          DO i=1,sNx
                0283           cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+sigma*cg2d_s(i,j,bi,bj)
                0284           cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-sigma*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0285          ENDDO
                0286         ENDDO
                0287        ENDDO
                0288       ENDDO
c90853988b Jean*0289 c     actualIts = 1
23bce0bbb8 Mart*0290 
                0291       CALL EXCH_S3D_RL( cg2d_r,1, myThid )
                0292 
                0293 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c90853988b Jean*0294       DO 10 it2d=1, numIters-1
                0295 c      actualIts = it2d
23bce0bbb8 Mart*0296 
                0297 C--    Solve preconditioning equation and update
                0298 C--    conjugate direction vector "s".
                0299 C--    z = M^-1 r
                0300        DO bj=myByLo(myThid),myByHi(myThid)
                0301         DO bi=myBxLo(myThid),myBxHi(myThid)
56f3d2cc76 Mart*0302 #ifdef TARGET_NEC_SX
                0303 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0304 #endif /* TARGET_NEC_SX */
c90853988b Jean*0305          DO j=1,sNy
                0306           DO i=1,sNx
                0307            cg2d_y(i,j,bi,bj) =
                0308      &       pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0309      &      +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0310      &      +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0311      &      +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0312      &      +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
23bce0bbb8 Mart*0313           ENDDO
                0314          ENDDO
                0315         ENDDO
                0316        ENDDO
                0317 
                0318        CALL EXCH_S3D_RL( cg2d_y, 1, myThid )
ff02675122 Jean*0319 
23bce0bbb8 Mart*0320 C==    v = A.z
                0321 C--    eta_qr = <z,r>
                0322 C--    delta = <z,v>
                0323 C--    Do the error calcuation here to consolidate global reductions
                0324        DO bj=myByLo(myThid),myByHi(myThid)
                0325         DO bi=myBxLo(myThid),myBxHi(myThid)
                0326          eta_qrNtile(bi,bj) = 0. _d 0
                0327          deltaTile(bi,bj)   = 0. _d 0
                0328          errTile(bi,bj) = 0. _d 0
56f3d2cc76 Mart*0329 #ifdef TARGET_NEC_SX
                0330 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0331 #endif /* TARGET_NEC_SX */
c90853988b Jean*0332          DO j=1,sNy
                0333           DO i=1,sNx
                0334            cg2d_v(i,j,bi,bj) =
                0335      &       aW2d(i  ,j  ,bi,bj)*cg2d_y(i-1,j  ,bi,bj)
                0336      &      +aW2d(i+1,j  ,bi,bj)*cg2d_y(i+1,j  ,bi,bj)
                0337      &      +aS2d(i  ,j  ,bi,bj)*cg2d_y(i  ,j-1,bi,bj)
                0338      &      +aS2d(i  ,j+1,bi,bj)*cg2d_y(i  ,j+1,bi,bj)
                0339      &      +aC2d(i  ,j  ,bi,bj)*cg2d_y(i  ,j  ,bi,bj)
23bce0bbb8 Mart*0340            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
c90853988b Jean*0341      &      +cg2d_y(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0342            deltaTile(bi,bj) = deltaTile(bi,bj)
c90853988b Jean*0343      &      +cg2d_y(i,j,bi,bj)*cg2d_v(i,j,bi,bj)
23bce0bbb8 Mart*0344            errTile(bi,bj) = errTile(bi,bj)
c90853988b Jean*0345      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
23bce0bbb8 Mart*0346           ENDDO
                0347          ENDDO
5237154b93 Jean*0348          tiledVec(bi,bj,1) = eta_qrNtile(bi,bj)
                0349          tiledVec(bi,bj,2) = deltaTile(bi,bj)
                0350          tiledVec(bi,bj,3) = errTile(bi,bj)
23bce0bbb8 Mart*0351         ENDDO
ff02675122 Jean*0352        ENDDO
                0353 
5237154b93 Jean*0354        CALL GLOBAL_SUM_VECTOR_RL( 3, tiledVec, sumVec, myThid )
ff02675122 Jean*0355 
5237154b93 Jean*0356        eta_qrN = sumVec(1)
                0357        delta   = sumVec(2)
                0358        err_sq  = sumVec(3)
ff02675122 Jean*0359 
764aea3440 Jean*0360        IF ( printResidual ) THEN
                0361         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0362          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
c90853988b Jean*0363      &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
764aea3440 Jean*0364          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0365      &                       SQUEEZE_RIGHT, myThid )
                0366         ENDIF
23bce0bbb8 Mart*0367        ENDIF
c90853988b Jean*0368        IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
                0369        IF ( err_sq .LT. minResidualSq ) THEN
                0370 C-     Store lowest residual solution
                0371          minResidualSq = err_sq
                0372          nIterMin = it2d
                0373          DO bj=myByLo(myThid),myByHi(myThid)
                0374           DO bi=myBxLo(myThid),myBxHi(myThid)
                0375            DO j=1,sNy
                0376             DO i=1,sNx
                0377              cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0378             ENDDO
                0379            ENDDO
                0380           ENDDO
                0381          ENDDO
                0382        ENDIF
ff02675122 Jean*0383 
23bce0bbb8 Mart*0384        cgBeta   = eta_qrN/eta_qrNM1
                0385        eta_qrNM1 = eta_qrN
c90853988b Jean*0386        alpha = delta - (cgBeta*cgBeta)*alpha
23bce0bbb8 Mart*0387        sigma = eta_qrN/alpha
ff02675122 Jean*0388 
c90853988b Jean*0389 C==    Update simultaneously solution and residual vectors (and Iter number)
23bce0bbb8 Mart*0390        DO bj=myByLo(myThid),myByHi(myThid)
                0391         DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0392          DO j=1,sNy
                0393           DO i=1,sNx
                0394            cg2d_s(i,j,bi,bj) = cg2d_y(i,j,bi,bj)
                0395      &                       + cgBeta*cg2d_s(i,j,bi,bj)
                0396            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0397      &                       + sigma*cg2d_s(i,j,bi,bj)
                0398            cg2d_q(i,j,bi,bj) = cg2d_v(i,j,bi,bj)
                0399      &                       + cgBeta*cg2d_q(i,j,bi,bj)
                0400            cg2d_r(i,j,bi,bj) = cg2d_r(i,j,bi,bj)
                0401      &                       - sigma*cg2d_q(i,j,bi,bj)
23bce0bbb8 Mart*0402           ENDDO
                0403          ENDDO
                0404         ENDDO
                0405        ENDDO
c90853988b Jean*0406 c      actualIts = it2d + 1
23bce0bbb8 Mart*0407 
                0408        CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
                0409 
                0410    10 CONTINUE
c90853988b Jean*0411       DO bj=myByLo(myThid),myByHi(myThid)
                0412        DO bi=myBxLo(myThid),myBxHi(myThid)
                0413         errTile(bi,bj) = 0. _d 0
                0414         DO j=1,sNy
                0415          DO i=1,sNx
                0416           errTile(bi,bj) = errTile(bi,bj)
                0417      &                   + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0418          ENDDO
                0419         ENDDO
                0420        ENDDO
                0421       ENDDO
                0422       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
23bce0bbb8 Mart*0423    11 CONTINUE
                0424 
c90853988b Jean*0425       IF ( nIterMin.GE.0 .AND. err_sq .GT. minResidualSq ) THEN
                0426 C-    use the lowest residual solution (instead of current one = last residual)
                0427         DO bj=myByLo(myThid),myByHi(myThid)
                0428          DO bi=myBxLo(myThid),myBxHi(myThid)
                0429           DO j=1,sNy
                0430            DO i=1,sNx
                0431              cg2d_x(i,j,bi,bj) = cg2d_min(i,j,bi,bj)
                0432            ENDDO
                0433           ENDDO
                0434          ENDDO
                0435         ENDDO
                0436       ENDIF
                0437 
23bce0bbb8 Mart*0438       IF (cg2dNormaliseRHS) THEN
                0439 C--   Un-normalise the answer
                0440         DO bj=myByLo(myThid),myByHi(myThid)
                0441          DO bi=myBxLo(myThid),myBxHi(myThid)
c90853988b Jean*0442           DO j=1,sNy
                0443            DO i=1,sNx
                0444             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
23bce0bbb8 Mart*0445            ENDDO
                0446           ENDDO
                0447          ENDDO
                0448         ENDDO
                0449       ENDIF
                0450 
                0451 C--   Return parameters to caller
c90853988b Jean*0452       lastResidual = SQRT(err_sq)
                0453       numIters = it2d
                0454 c     numIters = actualIts
23bce0bbb8 Mart*0455 
                0456 #endif /* ALLOW_SRCG */
                0457       RETURN
                0458       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!