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

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

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

revision 1.34.6.4 by heimbach, Tue Jul 8 15:18:29 2003 UTC revision 1.55 by jmc, Fri May 11 23:28:10 2012 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef TARGET_NEC_SX
6    C     set a sensible default for the outer loop unrolling parameter that can
7    C     be overriden in the Makefile with the DEFINES macro or in CPP_OPTIONS.h
8    #ifndef CG2D_OUTERLOOPITERS
9    # define CG2D_OUTERLOOPITERS 10
10    #endif
11    #endif /* TARGET_NEC_SX */
12    
13  CBOP  CBOP
14  C     !ROUTINE: CG2D  C     !ROUTINE: CG2D
15  C     !INTERFACE:  C     !INTERFACE:
16        SUBROUTINE CG2D(          SUBROUTINE CG2D(
17       I                cg2d_b,       U                cg2d_b, cg2d_x,
18       U                cg2d_x,       O                firstResidual, minResidualSq, lastResidual,
19       O                firstResidual,       U                numIters, nIterMin,
      O                lastResidual,  
      U                numIters,  
20       I                myThid )       I                myThid )
21  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
22  C     *==========================================================*  C     *==========================================================*
23  C     | SUBROUTINE CG2D                                            C     | SUBROUTINE CG2D
24  C     | o Two-dimensional grid problem conjugate-gradient          C     | o Two-dimensional grid problem conjugate-gradient
25  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
26  C     *==========================================================*  C     *==========================================================*
27  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
28  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
29  C     | This implementation assumes A is a five-diagonal            C     | This implementation assumes A is a five-diagonal
30  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
31  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
32  C     | two-dimensional space.                                      C     | two-dimensional space.
33  C     | Notes:                                                      C     | Notes:
34  C     | ======                                                      C     | ======
35  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
36  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
37  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
38  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
39  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
40  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
41  C     | the threads to solve several different problems              C     | the threads to solve several different problems
42  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
43  C     *==========================================================*  C     *==========================================================*
44  C     \ev  C     \ev
45    
# Line 44  C     === Global data === Line 49  C     === Global data ===
49  #include "SIZE.h"  #include "SIZE.h"
50  #include "EEPARAMS.h"  #include "EEPARAMS.h"
51  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
52  #include "CG2D.h"  #include "CG2D.h"
 #include "SURFACE.h"  
53    
54  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
55  C     === Routine arguments ===  C     === Routine arguments ===
56  C     myThid    - Thread on which I am working.  C     cg2d_b    :: The source term or "right hand side" (output: normalised RHS)
57  C     cg2d_b    - The source term or "right hand side"  C     cg2d_x    :: The solution (input: first guess)
58  C     cg2d_x    - The solution  C     firstResidual :: the initial residual before any iterations
59  C     firstResidual - the initial residual before any iterations  C     minResidualSq :: the lowest residual reached (squared)
60  C     lastResidual  - the actual residual reached  C     lastResidual  :: the actual residual reached
61  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  :: Inp: the maximum number of iterations allowed
62  C                 Exit:  the actual number of iterations used  C                  Out: the actual number of iterations used
63    C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
64    C                  Out: iteration number corresponding to lowest residual
65    C     myThid    :: Thread on which I am working.
66        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68        _RL  firstResidual        _RL  firstResidual
69          _RL  minResidualSq
70        _RL  lastResidual        _RL  lastResidual
71        INTEGER numIters        INTEGER numIters
72          INTEGER nIterMin
73        INTEGER myThid        INTEGER myThid
74    
75  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
76  C     === Local variables ====  C     === Local variables ====
77  C     actualIts      - Number of iterations taken  C     bi, bj     :: tile index in X and Y.
78  C     actualResidual - residual  C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
79  C     bi          - Block index in X and Y.  C     actualIts  :: actual CG iteration number
80  C     bj  C     err_sq     :: Measure of the square of the residual of Ax - b.
81  C     eta_qrN     - Used in computing search directions  C     eta_qrN    :: Used in computing search directions; suffix N and NM1
82  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     denote current and previous iterations respectively.
83  C     cgBeta        previous iterations respectively.  C     cgBeta     :: coeff used to update conjugate direction vector "s".
84  C     alpha    C     alpha      :: coeff used to update solution & residual
85  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
86  C                   useful debuggin/trouble shooting diagnostic.  C                   debugging/trouble shooting diagnostic. For neumann problems
87  C                   For neumann problems sumRHS needs to be ~0.  C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
88  C                   or they converge at a non-zero residual.  C     cg2d_min   :: used to store solution corresponding to lowest residual.
89  C     err         - Measure of residual of Ax - b, usually the norm.  C     msgBuf     :: Informational/error message buffer
90  C     I, J, N     - Loop counters ( N counts CG iterations )        INTEGER bi, bj
91          INTEGER i, j, it2d
92        INTEGER actualIts        INTEGER actualIts
93        _RL    actualResidual        _RL    cg2dTolerance_sq
94        INTEGER bi, bj                      _RL    err_sq,  errTile(nSx,nSy)
95        INTEGER I, J, it2d        _RL    eta_qrN, eta_qrNtile(nSx,nSy)
       _RL    err  
       _RL    eta_qrN  
96        _RL    eta_qrNM1        _RL    eta_qrNM1
97        _RL    cgBeta        _RL    cgBeta
98        _RL    alpha        _RL    alpha,   alphaTile(nSx,nSy)
99        _RL    sumRHS        _RL    sumRHS,  sumRHStile(nSx,nSy)
100        _RL    rhsMax        _RL    rhsMax
101        _RL    rhsNorm        _RL    rhsNorm
102          _RL    cg2d_min(1:sNx,1:sNy,nSx,nSy)
103        INTEGER OLw  #ifdef CG2D_SINGLECPU_SUM
104        INTEGER OLe        _RL    localBuf(1:sNx,1:sNy,nSx,nSy)
105        INTEGER OLn  #endif
106        INTEGER OLs        CHARACTER*(MAX_LEN_MBUF) msgBuf
107        INTEGER exchWidthX        LOGICAL printResidual
       INTEGER exchWidthY  
       INTEGER myNz  
108  CEOP  CEOP
109    
110    C--   Initialise auxiliary constant, some output variable and inverter
111  CcnhDebugStarts        cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
112  C     CHARACTER*(MAX_LEN_FNAM) suff        minResidualSq = -1. _d 0
113  CcnhDebugEnds        eta_qrNM1     =  1. _d 0
   
   
 C--   Initialise inverter  
       eta_qrNM1 = 1. _d 0  
   
 CcnhDebugStarts  
 C     _EXCH_XY_R8( cg2d_b, myThid )  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )  
 C     suff = 'unnormalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 C     STOP  
 CcnhDebugEnds  
114    
115  C--   Normalise RHS  C--   Normalise RHS
116        rhsMax = 0. _d 0        rhsMax = 0. _d 0
117        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
118         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
119          DO J=1,sNy          DO j=1,sNy
120           DO I=1,sNx           DO i=1,sNx
121            cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*cg2dNorm            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
122            rhsMax = MAX(ABS(cg2d_b(I,J,bi,bj)),rhsMax)            rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
123           ENDDO           ENDDO
124          ENDDO          ENDDO
125         ENDDO         ENDDO
# Line 134  C--   Normalise RHS Line 127  C--   Normalise RHS
127    
128        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
129  C-  Normalise RHS :  C-  Normalise RHS :
130  #ifdef LETS_MAKE_JAM        _GLOBAL_MAX_RL( rhsMax, myThid )
 C     _GLOBAL_MAX_R8( rhsMax, myThid )  
       rhsMax=1.  
 #else  
       _GLOBAL_MAX_R8( rhsMax, myThid )  
 Catm  rhsMax=1.  
 #endif  
131        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
132        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
133        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
134         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
135          DO J=1,sNy          DO j=1,sNy
136           DO I=1,sNx           DO i=1,sNx
137            cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj)*rhsNorm            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
138            cg2d_x(I,J,bi,bj) = cg2d_x(I,J,bi,bj)*rhsNorm            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
139           ENDDO           ENDDO
140          ENDDO          ENDDO
141         ENDDO         ENDDO
# Line 157  C- end Normalise RHS Line 144  C- end Normalise RHS
144        ENDIF        ENDIF
145    
146  C--   Update overlaps  C--   Update overlaps
147        _EXCH_XY_R8( cg2d_b, myThid )        CALL EXCH_XY_RL( cg2d_x, myThid )
       _EXCH_XY_R8( cg2d_x, myThid )  
 CcnhDebugStarts  
 C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  
 C     suff = 'normalised'  
 C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)  
 CcnhDebugEnds  
148    
149  C--   Initial residual calculation  C--   Initial residual calculation
       err    = 0. _d 0  
       sumRHS = 0. _d 0  
150        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
151         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
152          DO J=1,sNy          IF ( nIterMin.GE.0 ) THEN
153           DO I=1,sNx           DO j=1,sNy
154            cg2d_s(I,J,bi,bj) = 0.            DO i=1,sNx
155            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -              cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
156       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)            ENDDO
157       &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)           ENDDO
158       &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)          ENDIF
159       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)          DO j=0,sNy+1
160       &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)           DO i=0,sNx+1
161       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)            cg2d_s(i,j,bi,bj) = 0.
162       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)           ENDDO
163       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)          ENDDO
164       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*          sumRHStile(bi,bj) = 0. _d 0
165       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm          errTile(bi,bj)    = 0. _d 0
166    #ifdef TARGET_NEC_SX
167    !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
168    #endif /* TARGET_NEC_SX */
169            DO j=1,sNy
170             DO i=1,sNx
171              cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
172         &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
173         &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
174         &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
175         &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
176         &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
177       &    )       &    )
178            err            = err            +  #ifdef CG2D_SINGLECPU_SUM
179       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)            localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
180            sumRHS            = sumRHS            +  #else
181       &     cg2d_b(I,J,bi,bj)            errTile(bi,bj)    = errTile(bi,bj)
182         &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
183              sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
184    #endif
185           ENDDO           ENDDO
186          ENDDO          ENDDO
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189  C     _EXCH_XY_R8( cg2d_r, myThid )        CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
190  #ifdef LETS_MAKE_JAM  #ifdef CG2D_SINGLECPU_SUM
191        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
192  #else        CALL GLOBAL_SUM_SINGLECPU_RL(cg2d_b, sumRHS, OLx, OLy, myThid)
       CALL EXCH_XY_RL( cg2d_r, myThid )  
 #endif  
 C     _EXCH_XY_R8( cg2d_s, myThid )  
 #ifdef LETS_MAKE_JAM  
       CALL EXCH_XY_O1_R8_JAM( cg2d_s )  
193  #else  #else
194        CALL EXCH_XY_RL( cg2d_s, myThid )        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
195          CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
196  #endif  #endif
197         _GLOBAL_SUM_R8( sumRHS, myThid )        actualIts = 0
198         _GLOBAL_SUM_R8( err   , myThid )        firstResidual = SQRT(err_sq)
199         err = SQRT(err)        IF ( nIterMin.GE.0 ) THEN
200         actualIts      = 0          nIterMin = 0
201         actualResidual = err          minResidualSq = err_sq
202                ENDIF
203         IF ( debugLevel .GE. debLevA ) THEN  
204          printResidual = .FALSE.
205          IF ( debugLevel .GE. debLevZero ) THEN
206          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
207          write(*,'(A,1P2E22.14)')' cg2d: Sum(rhs),rhsMax = ',          printResidual = printResidualFreq.GE.1
208       &                                  sumRHS,rhsMax          WRITE(standardmessageunit,'(A,1P2E22.14)')
209          _END_MASTER( )       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
210         ENDIF          _END_MASTER( myThid )
211  C     _BARRIER        ENDIF
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( )  
       firstResidual=actualResidual  
212    
213        IF ( err .LT. cg2dTolerance ) GOTO 11        IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
214    
215  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
216        DO 10 it2d=1, numIters        DO 10 it2d=1, numIters
217    
 CcnhDebugStarts  
 C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',  
 C    &  actualResidual  
 CcnhDebugEnds  
218  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
219  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
        eta_qrN = 0. _d 0  
220         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
221          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
222           DO J=1,sNy           eta_qrNtile(bi,bj) = 0. _d 0
223            DO I=1,sNx  #ifdef TARGET_NEC_SX
224             cg2d_q(I,J,bi,bj) =  !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
225       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)  #endif /* TARGET_NEC_SX */
226       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)           DO j=1,sNy
227       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)            DO i=1,sNx
228       &     +pS(I  ,J  ,bi,bj)*cg2d_r(I  ,J-1,bi,bj)             cg2d_q(i,j,bi,bj) =
229       &     +pS(I  ,J+1,bi,bj)*cg2d_r(I  ,J+1,bi,bj)       &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
230         &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
231         &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
232         &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
233         &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
234  CcnhDebugStarts  CcnhDebugStarts
235  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)  c          cg2d_q(i,j,bi,bj) = cg2d_r(j  ,j  ,bi,bj)
236  CcnhDebugEnds  CcnhDebugEnds
237             eta_qrN = eta_qrN  #ifdef CG2D_SINGLECPU_SUM
238       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)            localBuf(i,j,bi,bj) =
239         &      cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
240    #else
241               eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
242         &     +cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
243    #endif
244            ENDDO            ENDDO
245           ENDDO           ENDDO
246          ENDDO          ENDDO
247         ENDDO         ENDDO
248    
249         _GLOBAL_SUM_R8(eta_qrN, myThid)  #ifdef CG2D_SINGLECPU_SUM
250  CcnhDebugStarts         CALL GLOBAL_SUM_SINGLECPU_RL( localBuf,eta_qrN,0,0,myThid )
251  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' eta_qrN = ',eta_qrN  #else
252  CcnhDebugEnds         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
253    #endif
254         cgBeta   = eta_qrN/eta_qrNM1         cgBeta   = eta_qrN/eta_qrNM1
255  CcnhDebugStarts  CcnhDebugStarts
256  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta  c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
257    c    &            ' eta_qrN=', eta_qrN, ' beta=', cgBeta
258  CcnhDebugEnds  CcnhDebugEnds
259         eta_qrNM1 = eta_qrN         eta_qrNM1 = eta_qrN
260    
261         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
262          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
263           DO J=1,sNy           DO j=1,sNy
264            DO I=1,sNx            DO i=1,sNx
265             cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj)             cg2d_s(i,j,bi,bj) = cg2d_q(i,j,bi,bj)
266       &                       + cgBeta*cg2d_s(I,J,bi,bj)       &                       + cgBeta*cg2d_s(i,j,bi,bj)
267            ENDDO            ENDDO
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270         ENDDO         ENDDO
271    
272  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between processes.
273  C--    processes.         CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
 C      _EXCH_XY_R8( cg2d_s, myThid )  
 #ifdef LETS_MAKE_JAM  
       CALL EXCH_XY_O1_R8_JAM( cg2d_s )  
 #else  
       CALL EXCH_XY_RL( cg2d_s, myThid )  
 #endif  
274    
275  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
276  C==    q = A.s  C==    q = A.s
        alpha = 0. _d 0  
277         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
278          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
279           DO J=1,sNy           alphaTile(bi,bj) = 0. _d 0
280            DO I=1,sNx  #ifdef TARGET_NEC_SX
281             cg2d_q(I,J,bi,bj) =  !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
282       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)  #endif /* TARGET_NEC_SX */
283       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)           DO j=1,sNy
284       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)            DO i=1,sNx
285       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)             cg2d_q(i,j,bi,bj) =
286       &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
287       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
288       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
289       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
290       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*       &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
291       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  #ifdef CG2D_SINGLECPU_SUM
292            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            localBuf(i,j,bi,bj) = cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
293    #else
294              alphaTile(bi,bj) = alphaTile(bi,bj)
295         &                     + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
296    #endif
297            ENDDO            ENDDO
298           ENDDO           ENDDO
299          ENDDO          ENDDO
300         ENDDO         ENDDO
301         _GLOBAL_SUM_R8(alpha,myThid)  #ifdef CG2D_SINGLECPU_SUM
302           CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, alpha, 0, 0, myThid)
303    #else
304           CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
305    #endif
306  CcnhDebugStarts  CcnhDebugStarts
307  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha  c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
308    c    &            ' SUM(s*q)=', alpha, ' alpha=', eta_qrN/alpha
309  CcnhDebugEnds  CcnhDebugEnds
310         alpha = eta_qrN/alpha         alpha = eta_qrN/alpha
311  CcnhDebugStarts  
312  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C==    Update simultaneously solution and residual vectors (and Iter number)
 CcnhDebugEnds  
       
 C==    Update solution and residual vectors  
313  C      Now compute "interior" points.  C      Now compute "interior" points.
        err = 0. _d 0  
314         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
315          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
316           DO J=1,sNy           errTile(bi,bj) = 0. _d 0
317            DO I=1,sNx           DO j=1,sNy
318             cg2d_x(I,J,bi,bj)=cg2d_x(I,J,bi,bj)+alpha*cg2d_s(I,J,bi,bj)            DO i=1,sNx
319             cg2d_r(I,J,bi,bj)=cg2d_r(I,J,bi,bj)-alpha*cg2d_q(I,J,bi,bj)             cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
320             err = err+cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)             cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
321    #ifdef CG2D_SINGLECPU_SUM
322               localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
323    #else
324               errTile(bi,bj) = errTile(bi,bj)
325         &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
326    #endif
327            ENDDO            ENDDO
328           ENDDO           ENDDO
329          ENDDO          ENDDO
330         ENDDO         ENDDO
331           actualIts = it2d
332    
333         _GLOBAL_SUM_R8( err   , myThid )  #ifdef CG2D_SINGLECPU_SUM
334         err = SQRT(err)         CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
        actualIts      = it2d  
        actualResidual = err  
        IF ( err .LT. cg2dTolerance ) GOTO 11  
 C      _EXCH_XY_R8(cg2d_r, myThid )  
 #ifdef LETS_MAKE_JAM  
       CALL EXCH_XY_O1_R8_JAM( cg2d_r )  
335  #else  #else
336        CALL EXCH_XY_RL( cg2d_r, myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
337  #endif  #endif
338           IF ( printResidual ) THEN
339            IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
340             WRITE(msgBuf,'(A,I6,A,1PE21.14)')
341         &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
342             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
343         &                       SQUEEZE_RIGHT, myThid )
344            ENDIF
345           ENDIF
346           IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
347           IF ( err_sq .LT. minResidualSq ) THEN
348    C-     Store lowest residual solution
349             minResidualSq = err_sq
350             nIterMin = it2d
351             DO bj=myByLo(myThid),myByHi(myThid)
352              DO bi=myBxLo(myThid),myBxHi(myThid)
353               DO j=1,sNy
354                DO i=1,sNx
355                 cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
356                ENDDO
357               ENDDO
358              ENDDO
359             ENDDO
360           ENDIF
361    
362           CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
363    
364     10 CONTINUE     10 CONTINUE
365     11 CONTINUE     11 CONTINUE
366    
367          IF ( nIterMin.GE.0 .AND. err_sq .GT. minResidualSq ) THEN
368    C-    use the lowest residual solution (instead of current one = last residual)
369            DO bj=myByLo(myThid),myByHi(myThid)
370             DO bi=myBxLo(myThid),myBxHi(myThid)
371              DO j=1,sNy
372               DO i=1,sNx
373                 cg2d_x(i,j,bi,bj) = cg2d_min(i,j,bi,bj)
374               ENDDO
375              ENDDO
376             ENDDO
377            ENDDO
378          ENDIF
379    
380        IF (cg2dNormaliseRHS) THEN        IF (cg2dNormaliseRHS) THEN
381  C--   Un-normalise the answer  C--   Un-normalise the answer
382          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
383           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
384            DO J=1,sNy            DO j=1,sNy
385             DO I=1,sNx             DO i=1,sNx
386              cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm              cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
387             ENDDO             ENDDO
388            ENDDO            ENDDO
389           ENDDO           ENDDO
390          ENDDO          ENDDO
391        ENDIF        ENDIF
392    
 C     The following exchange was moved up to solve_for_pressure  
 C     for compatibility with TAMC.  
 C     _EXCH_XY_R8(cg2d_x, myThid )  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  
 c    & actualIts, actualResidual  
 c     _END_MASTER( )  
   
393  C--   Return parameters to caller  C--   Return parameters to caller
394        lastResidual=actualResidual        lastResidual = SQRT(err_sq)
395        numIters=actualIts        numIters = actualIts
396    
397  CcnhDebugStarts  CcnhDebugStarts
398  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )  c     _EXCH_XY_RL(cg2d_x, myThid )
399  C     err    = 0. _d 0  c     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
400  C     DO bj=myByLo(myThid),myByHi(myThid)  c     err_sq = 0. _d 0
401  C      DO bi=myBxLo(myThid),myBxHi(myThid)  c     DO bj=myByLo(myThid),myByHi(myThid)
402  C       DO J=1,sNy  c      DO bi=myBxLo(myThid),myBxHi(myThid)
403  C        DO I=1,sNx  c       DO j=1,sNy
404  C         cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -  c        DO i=1,sNx
405  C    &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)  c         cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
406  C    &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)  c    &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
407  C    &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)  c    &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
408  C    &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)  c    &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
409  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  c    &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
410  C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  c    &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
411  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  c    &    )
412  C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))  c         err_sq = err_sq + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
413  C         err            = err            +  c        ENDDO
414  C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  c       ENDDO
415  C        ENDDO  c      ENDDO
416  C       ENDDO  c     ENDDO
417  C      ENDDO  c     _GLOBAL_SUM_RL( err_sq, myThid )
418  C     ENDDO  c     write(*,*) 'cg2d: Ax - b = ',SQRT(err_sq)
 C     _GLOBAL_SUM_R8( err   , myThid )  
 C     write(*,*) 'cg2d: Ax - b = ',SQRT(err)  
419  CcnhDebugEnds  CcnhDebugEnds
420    
421        RETURN        RETURN
Legend:
Removed from v.1.34.6.4  
changed lines
  Added in v.1.55
  ViewVC Help
Powered by ViewVC 1.1.22