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

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

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

revision 1.1 by utke, Wed Jan 11 06:31:20 2006 UTC revision 1.3 by utke, Thu Jan 11 15:19:59 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 50  C     === Global data === Line 50  C     === Global data ===
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 131  C     _GLOBAL_MAX_R8( rhsMax, myThid ) Line 131  C     _GLOBAL_MAX_R8( rhsMax, myThid )
131        rhsMax=1.        rhsMax=1.
132  #else  #else
133        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_R8( rhsMax, myThid )
 Catm  rhsMax=1.  
134  #endif  #endif
135        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
136        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 165  C--   Initial residual calculation
165          errTile    = 0. _d 0          errTile    = 0. _d 0
166          DO J=1,sNy          DO J=1,sNy
167           DO I=1,sNx           DO I=1,sNx
168    c         ks = ksurfC(I,J,bi,bj)
169            cg2d_s(I,J,bi,bj) = 0.            cg2d_s(I,J,bi,bj) = 0.
170            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -            cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -
171       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)       &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)
# 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 227  c     _END_MASTER( myThid ) Line 226  c     _END_MASTER( myThid )
226  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
227        it2d=1        it2d=1
228                
229        DO while ((it2d1 .le. numIters .and. err .ge. cg2dTolerance))        DO while ((it2d .le. numIters .and. err .ge. cg2dTolerance))
230    
231  CcnhDebugStarts  CcnhDebugStarts
232  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',
# 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)
# Line 321  CcnhDebugEnds Line 321  CcnhDebugEnds
321  CcnhDebugStarts  CcnhDebugStarts
322  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
323  CcnhDebugEnds  CcnhDebugEnds
324        
325  C==    Update solution and residual vectors  C==    Update solution and residual vectors
326  C      Now compute "interior" points.  C      Now compute "interior" points.
327         err = 0. _d 0         err = 0. _d 0
# Line 349  C      _EXCH_XY_R8(cg2d_r, myThid ) Line 349  C      _EXCH_XY_R8(cg2d_r, myThid )
349  #else  #else
350        CALL EXCH_XY_RL( cg2d_r, myThid )        CALL EXCH_XY_RL( cg2d_r, myThid )
351  #endif  #endif
352           it2d=it2d+1
353    
354        end do        end do
355        end if        end if
# Line 370  C     The following exchange was moved u Line 371  C     The following exchange was moved u
371  C     for compatibility with TAMC.  C     for compatibility with TAMC.
372  C     _EXCH_XY_R8(cg2d_x, myThid )  C     _EXCH_XY_R8(cg2d_x, myThid )
373  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
374  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
375  c    & actualIts, actualResidual  c    & actualIts, actualResidual
376  c     _END_MASTER( myThid )  c     _END_MASTER( myThid )
377    
# Line 394  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I Line 395  C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I
395  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)
396  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
397  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))
398  C         err            = err            +  C         err            = err            +
399  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)
400  C        ENDDO  C        ENDDO
401  C       ENDDO  C       ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22