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

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

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

revision 1.7 by jmc, Thu Mar 8 20:45:53 2001 UTC revision 1.15 by jmc, Fri Feb 4 19:30:33 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  #define VERBOSE  #define VERBOSE
8    
9    CBOP
10    C     !ROUTINE: CG3D
11    C     !INTERFACE:
12        SUBROUTINE CG3D(          SUBROUTINE CG3D(  
13         I                cg3d_b,
14         U                cg3d_x,
15         O                firstResidual,
16         O                lastResidual,
17         U                numIters,
18       I                myThid )       I                myThid )
19  C     /==========================================================\  C     !DESCRIPTION: \bv
20  C     | SUBROUTINE CG3D                                          |  C     *==========================================================*
21  C     | o Three-dimensional grid problem conjugate-gradient      |  C     | SUBROUTINE CG3D                                          
22  C     |   inverter (with preconditioner).                        |  C     | o Three-dimensional grid problem conjugate-gradient      
23  C     |==========================================================|  C     |   inverter (with preconditioner).                        
24  C     | Con. grad is an iterative procedure for solving Ax = b.  |  C     *==========================================================*
25  C     | It requires the A be symmetric.                          |  C     | Con. grad is an iterative procedure for solving Ax = b.  
26  C     | This implementation assumes A is a five-diagonal         |  C     | It requires the A be symmetric.                          
27  C     | matrix of the form that arises in the discrete           |  C     | This implementation assumes A is a seven-diagonal          
28  C     | representation of the del^2 operator in a                |  C     | matrix of the form that arises in the discrete            
29  C     | two-dimensional space.                                   |  C     | representation of the del^2 operator in a                
30  C     | Notes:                                                   |  C     | three-dimensional space.                                    
31  C     | ======                                                   |  C     | Notes:                                                    
32  C     | This implementation can support shared-memory            |  C     | ======                                                    
33  C     | multi-threaded execution. In order to do this COMMON     |  C     | This implementation can support shared-memory              
34  C     | blocks are used for many of the arrays - even ones that  |  C     | multi-threaded execution. In order to do this COMMON      
35  C     | are only used for intermedaite results. This design is   |  C     | blocks are used for many of the arrays - even ones that    
36  C     | OK if you want to all the threads to collaborate on      |  C     | are only used for intermedaite results. This design is    
37  C     | solving the same problem. On the other hand if you want  |  C     | OK if you want to all the threads to collaborate on        
38  C     | the threads to solve several different problems          |  C     | solving the same problem. On the other hand if you want    
39  C     | concurrently this implementation will not work.          |  C     | the threads to solve several different problems            
40  C     \==========================================================/  C     | concurrently this implementation will not work.          
41        IMPLICIT NONE  C     *==========================================================*
42    C     \ev
43    
44    C     !USES:
45          IMPLICIT NONE
46  C     === Global data ===  C     === Global data ===
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 38  C     === Global data === Line 50  C     === Global data ===
50  #include "GRID.h"  #include "GRID.h"
51  #include "CG3D.h"  #include "CG3D.h"
52    
53    C     !INPUT/OUTPUT PARAMETERS:
54  C     === Routine arguments ===  C     === Routine arguments ===
55  C     myThid - Thread on which I am working.  C     myThid    - Thread on which I am working.
56    C     cg2d_b    - The source term or "right hand side"
57    C     cg2d_x    - The solution
58    C     firstResidual - the initial residual before any iterations
59    C     lastResidual  - the actual residual reached
60    C     numIters  - Entry: the maximum number of iterations allowed
61    C                 Exit:  the actual number of iterations used
62          _RL  cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
63          _RL  cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
64          _RL  firstResidual
65          _RL  lastResidual
66          INTEGER numIters
67        INTEGER myThid        INTEGER myThid
68    
69    
70  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
71    
72    C     !LOCAL VARIABLES:
73  C     === Local variables ====  C     === Local variables ====
74  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
75  C     actualResidual - residual  C     actualResidual - residual
# Line 64  C     I, J, N     - Loop counters ( N co Line 90  C     I, J, N     - Loop counters ( N co
90        INTEGER bi, bj                      INTEGER bi, bj              
91        INTEGER I, J, K, it3d        INTEGER I, J, K, it3d
92        INTEGER KM1, KP1        INTEGER KM1, KP1
93        _RL    err        _RL    err, errTile
94        _RL    eta_qrN        _RL    eta_qrN, eta_qrNtile
95        _RL    eta_qrNM1        _RL    eta_qrNM1
96        _RL    cgBeta        _RL    cgBeta
97        _RL    alpha        _RL    alpha , alphaTile
98        _RL    sumRHS        _RL    sumRHS, sumRHStile
99        _RL    rhsMax        _RL    rhsMax
100        _RL    rhsNorm        _RL    rhsNorm
101    
# Line 81  C     I, J, N     - Loop counters ( N co Line 107  C     I, J, N     - Loop counters ( N co
107        INTEGER exchWidthY        INTEGER exchWidthY
108        INTEGER myNz        INTEGER myNz
109        _RL     topLevTerm        _RL     topLevTerm
110    CEOP
111    
112    ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
113    
114    
115  C--   Initialise inverter  C--   Initialise inverter
116        eta_qrNM1 = 1. D0        eta_qrNM1 = 1. D0
# Line 124  C--   Initial residual calculation (with Line 154  C--   Initial residual calculation (with
154        sumRHS = 0. _d 0        sumRHS = 0. _d 0
155        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
156         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
157            errTile    = 0. _d 0
158            sumRHStile = 0. _d 0
159          DO K=1,Nr          DO K=1,Nr
160           KM1 = K-1           KM1 = K-1
161           IF ( K .EQ. 1 ) KM1 = 1           IF ( K .EQ. 1 ) KM1 = 1
# Line 150  C--   Initial residual calculation (with Line 182  C--   Initial residual calculation (with
182       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
183       &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)       &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)
184       &     )       &     )
185             err = err             errTile = errTile
186       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
187             sumRHS = sumRHS             sumRHStile = sumRHStile
188       &     +cg3d_b(I,J,K,bi,bj)       &     +cg3d_b(I,J,K,bi,bj)
189            ENDDO            ENDDO
190           ENDDO           ENDDO
191          ENDDO          ENDDO
192            err    = err    + errTile
193            sumRHS = sumRHS + sumRHStile
194         ENDDO         ENDDO
195        ENDDO        ENDDO
196  C     _EXCH_XYZ_R8( cg3d_r, myThid )  C     _EXCH_XYZ_R8( cg3d_r, myThid )
# Line 186  C     _EXCH_XYZ_R8( cg3d_s, myThid ) Line 220  C     _EXCH_XYZ_R8( cg3d_s, myThid )
220        _GLOBAL_SUM_R8( sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHS, myThid )
221        _GLOBAL_SUM_R8( err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
222                
223        _BEGIN_MASTER( myThid )        IF ( debugLevel .GE. debLevZero ) THEN
224        write(0,'(A,1PE30.14)') ' cg3d: Sum(rhs) = ',sumRHS          _BEGIN_MASTER( myThid )
225        _END_MASTER( )          write(*,'(A,1P2E22.14)')
226         &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
227            _END_MASTER( myThid )
228          ENDIF
229    
230        actualIts      = 0        actualIts      = 0
231        actualResidual = SQRT(err)        actualResidual = SQRT(err)
232  C     _BARRIER  C     _BARRIER
233        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
234         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
235       &  actualIts, actualResidual  c    &  actualIts, actualResidual
236        _END_MASTER( )  c     _END_MASTER( myThid )
237          firstResidual=actualResidual
238    
239  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, cg3dMaxIters
241    
242  CcnhDebugStarts  CcnhDebugStarts
243  #ifdef VERBOSE  #ifdef VERBOSE
244         IF ( mod(it3d-1,10).EQ.0)  c      IF ( mod(it3d-1,10).EQ.0)
245       &   WRITE(0,*) ' CG3D: Iteration ',it3d-1,  c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,
246       &      ' residual = ',actualResidual  c    &      ' residual = ',actualResidual
247  #endif  #endif
248  CcnhDebugEnds  CcnhDebugEnds
249         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
# Line 218  C            want eta_qrN for the interi Line 256  C            want eta_qrN for the interi
256         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
257         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
258          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
259             eta_qrNtile = 0. _d 0
260           DO K=1,1           DO K=1,1
261            DO J=1-1,sNy+1            DO J=1-1,sNy+1
262             DO I=1-1,sNx+1             DO I=1-1,sNx+1
# Line 247  caja       ENDDO Line 286  caja       ENDDO
286  caja      ENDIF  caja      ENDIF
287            DO J=1,sNy            DO J=1,sNy
288             DO I=1,sNx             DO I=1,sNx
289              eta_qrN = eta_qrN              eta_qrNtile = eta_qrNtile
290       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
291             ENDDO             ENDDO
292            ENDDO            ENDDO
# Line 262  caja      ENDIF Line 301  caja      ENDIF
301            ENDDO            ENDDO
302            DO J=1,sNy            DO J=1,sNy
303             DO I=1,sNx             DO I=1,sNx
304              eta_qrN = eta_qrN              eta_qrNtile = eta_qrNtile
305       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
306             ENDDO             ENDDO
307            ENDDO            ENDDO
308           ENDDO           ENDDO
309             eta_qrN = eta_qrN + eta_qrNtile
310          ENDDO          ENDDO
311         ENDDO         ENDDO
312  caja  caja
# Line 287  caja Line 327  caja
327    
328         _GLOBAL_SUM_R8(eta_qrN, myThid)         _GLOBAL_SUM_R8(eta_qrN, myThid)
329  CcnhDebugStarts  CcnhDebugStarts
330  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
331  CcnhDebugEnds  CcnhDebugEnds
332         cgBeta   = eta_qrN/eta_qrNM1         cgBeta   = eta_qrN/eta_qrNM1
333  CcnhDebugStarts  CcnhDebugStarts
334  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta
335  CcnhDebugEnds  CcnhDebugEnds
336         eta_qrNM1 = eta_qrN         eta_qrNM1 = eta_qrN
337    
# Line 315  C==    q = A.s Line 355  C==    q = A.s
355       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom
356         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
357          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
358             alphaTile = 0. _d 0
359           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
360            DO K=1,1            DO K=1,1
361             DO J=1,sNy             DO J=1,sNy
# Line 331  C==    q = A.s Line 372  C==    q = A.s
372       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
373       &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
374       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
375               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alphaTile = alphaTile
376         &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
377              ENDDO              ENDDO
378             ENDDO             ENDDO
379            ENDDO            ENDDO
# Line 349  C==    q = A.s Line 391  C==    q = A.s
391       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
392       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
393       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
394               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alphaTile = alphaTile
395         &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
396              ENDDO              ENDDO
397             ENDDO             ENDDO
398            ENDDO            ENDDO
# Line 370  C==    q = A.s Line 413  C==    q = A.s
413       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
414       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
415       &     -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
416              alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)              alphaTile = alphaTile
417         &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
418             ENDDO             ENDDO
419            ENDDO            ENDDO
420           ENDDO           ENDDO
# Line 389  C==    q = A.s Line 433  C==    q = A.s
433       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
434       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
435       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
436               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alphaTile = alphaTile
437         &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
438              ENDDO              ENDDO
439             ENDDO             ENDDO
440            ENDDO            ENDDO
441           ENDIF           ENDIF
442             alpha = alpha + alphaTile
443          ENDDO          ENDDO
444         ENDDO         ENDDO
445         _GLOBAL_SUM_R8(alpha,myThid)         _GLOBAL_SUM_R8(alpha,myThid)
446  CcnhDebugStarts  CcnhDebugStarts
447  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
448  CcnhDebugEnds  CcnhDebugEnds
449         alpha = eta_qrN/alpha         alpha = eta_qrN/alpha
450  CcnhDebugStarts  CcnhDebugStarts
451  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
452  CcnhDebugEnds  CcnhDebugEnds
453            
454  C==    Update solution and residual vectors  C==    Update solution and residual vectors
# Line 410  C      Now compute "interior" points. Line 456  C      Now compute "interior" points.
456         err = 0. _d 0         err = 0. _d 0
457         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
458          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
459            errTile    = 0. _d 0
460           DO K=1,Nr           DO K=1,Nr
461            DO J=1,sNy            DO J=1,sNy
462             DO I=1,sNx             DO I=1,sNx
# Line 417  C      Now compute "interior" points. Line 464  C      Now compute "interior" points.
464       &            +alpha*cg3d_s(I,J,K,bi,bj)       &            +alpha*cg3d_s(I,J,K,bi,bj)
465              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)
466       &            -alpha*cg3d_q(I,J,K,bi,bj)       &            -alpha*cg3d_q(I,J,K,bi,bj)
467              err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)             errTile = errTile
468         &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
469             ENDDO             ENDDO
470            ENDDO            ENDDO
471           ENDDO           ENDDO
472             err = err + errTile
473          ENDDO          ENDDO
474         ENDDO         ENDDO
475    
# Line 459  C--   Un-normalise the answer Line 508  C--   Un-normalise the answer
508        ENDDO        ENDDO
509    
510  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )
511        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
512         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
513       &  actualIts, actualResidual  c    &  actualIts, actualResidual
514        _END_MASTER( )  c     _END_MASTER( myThid )
515          lastResidual=actualResidual
516          numIters=actualIts
517    
518  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
519    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22