/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_shallow_openad3/cg2d.F
ViewVC logotype

Diff of /MITgcm_contrib/heimbach/OpenAD/code_shallow_openad3/cg2d.F

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

revision 1.1 by utke, Tue Aug 1 15:44:28 2006 UTC revision 1.2 by utke, Thu May 10 01:51:48 2007 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: CG2D  C     !ROUTINE: CG2D
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE CG2D(          SUBROUTINE CG2D(
10       I                cg2d_b,       I                cg2d_b,
11       U                cg2d_x,       U                cg2d_x,
12       O                firstResidual,       O                firstResidual,
# Line 15  C     !INTERFACE: Line 15  C     !INTERFACE:
15       I                myThid )       I                myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CG2D                                            C     | SUBROUTINE CG2D
19  C     | o Two-dimensional grid problem conjugate-gradient          C     | o Two-dimensional grid problem conjugate-gradient
20  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
21  C     *==========================================================*  C     *==========================================================*
22  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
23  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
24  C     | This implementation assumes A is a five-diagonal            C     | This implementation assumes A is a five-diagonal
25  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
26  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
27  C     | two-dimensional space.                                      C     | two-dimensional space.
28  C     | Notes:                                                      C     | Notes:
29  C     | ======                                                      C     | ======
30  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
31  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
32  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
33  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
34  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
35  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
36  C     | the threads to solve several different problems              C     | the threads to solve several different problems
37  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
38  C     *==========================================================*  C     *==========================================================*
39  C     \ev  C     \ev
40    
# Line 44  C     === Global data === Line 44  C     === Global data ===
44  #include "SIZE.h"  #include "SIZE.h"
45  #include "EEPARAMS.h"  #include "EEPARAMS.h"
46  #include "PARAMS.h"  #include "PARAMS.h"
 #include "GRID.h"  
47  #include "CG2D.h"  #include "CG2D.h"
48  #include "SURFACE.h"  c#include "GRID.h"
49    c#include "SURFACE.h"
50    
51  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
52  C     === Routine arguments ===  C     === Routine arguments ===
53  C     myThid    - Thread on which I am working.  C     myThid    :: Thread on which I am working.
54  C     cg2d_b    - The source term or "right hand side"  C     cg2d_b    :: The source term or "right hand side"
55  C     cg2d_x    - The solution  C     cg2d_x    :: The solution
56  C     firstResidual - the initial residual before any iterations  C     firstResidual :: the initial residual before any iterations
57  C     lastResidual  - the actual residual reached  C     lastResidual  :: the actual residual reached
58  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  :: Entry: the maximum number of iterations allowed
59  C                 Exit:  the actual number of iterations used  C                  Exit:  the actual number of iterations used
60        _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)
61        _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)
62        _RL  firstResidual        _RL  firstResidual
# Line 66  C                 Exit:  the actual numb Line 66  C                 Exit:  the actual numb
66    
67  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
68  C     === Local variables ====  C     === Local variables ====
69  C     actualIts      - Number of iterations taken  C     actualIts      :: Number of iterations taken
70  C     actualResidual - residual  C     actualResidual :: residual
71  C     bi          - Block index in X and Y.  C     bi, bj     :: Block index in X and Y.
72  C     bj  C     eta_qrN    :: Used in computing search directions
 C     eta_qrN     - Used in computing search directions  
73  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
74  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
75  C     alpha    C     alpha
76  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a
77  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
78  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
79  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
80  C     err         - Measure of residual of Ax - b, usually the norm.  C     err        :: Measure of residual of Ax - b, usually the norm.
81  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, it2d :: Loop counters ( it2d counts CG iterations )
82        INTEGER actualIts        INTEGER actualIts
83        _RL    actualResidual        _RL    actualResidual
84        INTEGER bi, bj                      INTEGER bi, bj
85        INTEGER I, J, it2d        INTEGER I, J, it2d
86    c     INTEGER ks
87        _RL    err,errTile        _RL    err,errTile
88        _RL    eta_qrN,eta_qrNtile        _RL    eta_qrN,eta_qrNtile
89        _RL    eta_qrNM1        _RL    eta_qrNM1
# Line 92  C     I, J, N     - Loop counters ( N co Line 92  C     I, J, N     - Loop counters ( N co
92        _RL    sumRHS,sumRHStile        _RL    sumRHS,sumRHStile
93        _RL    rhsMax        _RL    rhsMax
94        _RL    rhsNorm        _RL    rhsNorm
95          CHARACTER*(MAX_LEN_MBUF) msgBuf
96  CEOP  CEOP
97    
98    
# Line 131  C     _GLOBAL_MAX_R8( rhsMax, myThid ) Line 132  C     _GLOBAL_MAX_R8( rhsMax, myThid )
132        rhsMax=1.        rhsMax=1.
133  #else  #else
134        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_R8( rhsMax, myThid )
 Catm  rhsMax=1.  
135  #endif  #endif
136        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
137        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
# Line 166  C--   Initial residual calculation Line 166  C--   Initial residual calculation
166          errTile    = 0. _d 0          errTile    = 0. _d 0
167          DO J=1,sNy          DO J=1,sNy
168           DO I=1,sNx           DO I=1,sNx
169    c         ks = ksurfC(I,J,bi,bj)
170            cg2d_s(I,J,bi,bj) = 0.            cg2d_s(I,J,bi,bj) = 0.
171            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -            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)       &    (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)       &    +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)       &    +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)       &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)
176       &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    +aC2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
      &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
      &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  
      &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  
177       &    )       &    )
178            errTile        = errTile        +  c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
179       &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
180            sumRHStile        = sumRHStile        +  c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
181       &     cg2d_b(I,J,bi,bj)  c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
182    c    &    -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
183    c    &                        *cg2d_x(I  ,J  ,bi,bj)
184    c    &                        /deltaTMom/deltaTfreesurf*cg2dNorm
185    c    &    )
186              errTile    = errTile + cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
187              sumRHStile = sumRHStile + cg2d_b(I,J,bi,bj)
188           ENDDO           ENDDO
189          ENDDO          ENDDO
190          sumRHS = sumRHS + sumRHStile          sumRHS = sumRHS + sumRHStile
191          err    = err    + errTile          err    = err    + errTile
192         ENDDO         ENDDO
193        ENDDO        ENDDO
 C     _EXCH_XY_R8( cg2d_r, myThid )  
194  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
195        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
196  #else  #else
197        CALL EXCH_XY_RL( cg2d_r, myThid )        CALL EXCH_XY_RL( cg2d_r, myThid )
198  #endif  #endif
 C     _EXCH_XY_R8( cg2d_s, myThid )  
199  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
200        CALL EXCH_XY_O1_R8_JAM( cg2d_s )        CALL EXCH_XY_O1_R8_JAM( cg2d_s )
201  #else  #else
# Line 209  C     _EXCH_XY_R8( cg2d_s, myThid ) Line 209  C     _EXCH_XY_R8( cg2d_s, myThid )
209    
210         IF ( debugLevel .GE. debLevZero ) THEN         IF ( debugLevel .GE. debLevZero ) THEN
211          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
212          write(standardmessageunit,'(A,1P2E22.14)')          WRITE(standardmessageunit,'(A,1P2E22.14)')
213       &  ' cg2d: Sum(rhs),rhsMax = ',       &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
      &                                  sumRHS,rhsMax  
214          _END_MASTER( myThid )          _END_MASTER( myThid )
215         ENDIF         ENDIF
216  C     _BARRIER  C     _BARRIER
217  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
218  c      WRITE(standardmessageunit,'(A,I6,1PE30.14)')  c      WRITE(standardmessageunit,'(A,I6,1PE30.14)')
219  c    & ' CG2D iters, err = ',  c    & ' CG2D iters, err = ',
220  c    & actualIts, actualResidual  c    & actualIts, actualResidual
221  c     _END_MASTER( myThid )  c     _END_MASTER( myThid )
222        firstResidual=actualResidual        firstResidual=actualResidual
# Line 241  C--    conjugate direction vector "s". Line 240  C--    conjugate direction vector "s".
240           eta_qrNtile = 0. _d 0           eta_qrNtile = 0. _d 0
241           DO J=1,sNy           DO J=1,sNy
242            DO I=1,sNx            DO I=1,sNx
243             cg2d_q(I,J,bi,bj) =             cg2d_q(I,J,bi,bj) =
244       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)       &      pC(I  ,J  ,bi,bj)*cg2d_r(I  ,J  ,bi,bj)
245       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)       &     +pW(I  ,J  ,bi,bj)*cg2d_r(I-1,J  ,bi,bj)
246       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)       &     +pW(I+1,J  ,bi,bj)*cg2d_r(I+1,J  ,bi,bj)
# Line 296  C==    q = A.s Line 295  C==    q = A.s
295           alphaTile = 0. _d 0           alphaTile = 0. _d 0
296           DO J=1,sNy           DO J=1,sNy
297            DO I=1,sNx            DO I=1,sNx
298             cg2d_q(I,J,bi,bj) =  c          ks = ksurfC(I,J,bi,bj)
299               cg2d_q(I,J,bi,bj) =
300       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)       &     aW2d(I  ,J  ,bi,bj)*cg2d_s(I-1,J  ,bi,bj)
301       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)       &    +aW2d(I+1,J  ,bi,bj)*cg2d_s(I+1,J  ,bi,bj)
302       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)       &    +aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J-1,bi,bj)
303       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)       &    +aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J+1,bi,bj)
304       &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    +aC2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
305       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aW2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
306       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
307       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)  c    &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
308       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*  c    &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
309       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTfreesurf*cg2dNorm  c    &    -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)*recip_Bo(i,j,bi,bj)
310    c    &                        *cg2d_s(I  ,J  ,bi,bj)
311    c    &                        /deltaTMom/deltaTfreesurf*cg2dNorm
312            alphaTile = alphaTile+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            alphaTile = alphaTile+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
313            ENDDO            ENDDO
314           ENDDO           ENDDO
# Line 321  CcnhDebugEnds Line 323  CcnhDebugEnds
323  CcnhDebugStarts  CcnhDebugStarts
324  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
325  CcnhDebugEnds  CcnhDebugEnds
326        
327  C==    Update solution and residual vectors  C==    Update solution and residual vectors
328  C      Now compute "interior" points.  C      Now compute "interior" points.
329         err = 0. _d 0         err = 0. _d 0
# Line 343  C      Now compute "interior" points. Line 345  C      Now compute "interior" points.
345         err = SQRT(err)         err = SQRT(err)
346         actualIts      = it2d         actualIts      = it2d
347         actualResidual = err         actualResidual = err
348  C      _EXCH_XY_R8(cg2d_r, myThid )         IF ( debugLevel.GT.debLevB ) THEN
349    c       IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
350    c    &     ) THEN
351            _BEGIN_MASTER( myThid )
352            WRITE(msgBuf,'(A,I6,A,1PE21.14)')
353         &    ' cg2d: iter=', actualIts, ' ; resid.= ', actualResidual
354            CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
355            _END_MASTER( myThid )
356    c       ENDIF
357           ENDIF
358  #ifdef LETS_MAKE_JAM  #ifdef LETS_MAKE_JAM
359        CALL EXCH_XY_O1_R8_JAM( cg2d_r )        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
360  #else  #else
# Line 371  C     The following exchange was moved u Line 382  C     The following exchange was moved u
382  C     for compatibility with TAMC.  C     for compatibility with TAMC.
383  C     _EXCH_XY_R8(cg2d_x, myThid )  C     _EXCH_XY_R8(cg2d_x, myThid )
384  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
385  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
386  c    & actualIts, actualResidual  c    & actualIts, actualResidual
387  c     _END_MASTER( myThid )  c     _END_MASTER( myThid )
388    
# Line 395  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I Line 406  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I
406  C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
407  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
408  C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))  C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))
409  C         err            = err            +  C         err            = err            +
410  C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
411  C        ENDDO  C        ENDDO
412  C       ENDDO  C       ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22