Changeset 787


Ignore:
Timestamp:
2008-01-11T14:33:42+01:00 (13 years ago)
Author:
rblod
Message:

Improve PCG algortithm in MPI case, see ticket #47

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/SOL/solpcg.F90

    r719 r787  
    5252      !!      where q is the preconditioning matrix = diagonal matrix of the 
    5353      !!                                              diagonal elements of a 
    54       !!      Initialization: 
     54      !!      Initialization  : 
    5555      !!         x(o) = gcx 
    5656      !!         r(o) = d(o) = pgcb - pa.x(o) 
    5757      !!         rr(o)= < r(o) , r(o) >_q 
    58       !!      Iteration n   : 
    59       !!         z(n)   = pa.d(n) 
    60       !!         alp(n) = rr(n) / < z(n) , d(n) >_q 
     58      !!      Iteration 1     : 
     59      !!         standard PCG algorithm 
     60      !!      Iteration n > 1 : 
     61      !!         s(n)   = pa.r(n) 
     62      !!         gam(n) = < r(n) , r(n) >_q 
     63      !!         del(n) = < r(n) , s(n) >_q 
     64      !!         bet(n) = gam(n) / gam(n-1) 
     65      !!         d(n)   = r(n) + bet(n) d(n-1) 
     66      !!         z(n)   = s(n) + bet(n) z(n-1)  
     67      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1)  
     68      !!         alp(n) = gam(n) / sig(n)  
    6169      !!         x(n+1) = x(n) + alp(n) d(n) 
    6270      !!         r(n+1) = r(n) - alp(n) z(n) 
    63       !!         rr(n+1)= < r(n+1) , r(n+1) >_q 
    64       !!         bet(n) = rr(n+1) / rr(n) 
    65       !!         r(n+1) = r(n+1) + bet(n+1) d(n) 
    6671      !!      Convergence test : 
    6772      !!         rr(n+1) / < gcb , gcb >_q   =< epsr 
     
    7378      !! References : 
    7479      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6. 
     80      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U. 
    7581      !! 
    7682      !! History : 
     
    8288      !!        !  96-11  (A. Weaver)  correction to preconditioning 
    8389      !!   8.5  !  02-08  (G. Madec)  F90: Free form 
     90      !!        !  08-01  (R. Benshila) mpp optimization 
    8491      !!---------------------------------------------------------------------- 
    8592      !! * Arguments 
     
    9198      !! * Local declarations 
    9299      INTEGER ::   ji, jj, jn                ! dummy loop indices 
    93       REAL(wp) ::   zgcad                    ! temporary scalars 
     100      REAL(wp) ::  zgcad                     ! temporary scalars 
     101      REAL(wp), DIMENSION(2) :: zsum 
     102      REAL(wp), DIMENSION(jpi,jpj) :: zgcr 
    94103      !!---------------------------------------------------------------------- 
    95104 
     105      ! Initialization of the algorithm with standard PCG 
     106      ! ------------------------------------------------- 
     107 
     108      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition 
     109 
     110      ! gcr   = gcb-a.gcx 
     111      ! gcdes = gcr 
     112 
     113      DO jj = 2, jpjm1 
     114         DO ji = fs_2, fs_jpim1   ! vector opt. 
     115            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   & 
     116               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   & 
     117               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   & 
     118               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   & 
     119               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   ) 
     120            gcr  (ji,jj) = zgcad 
     121            gcdes(ji,jj) = zgcad 
     122         END DO 
     123      END DO 
     124 
     125      ! rnorme = (gcr,gcr) 
     126      rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  ) 
     127      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain 
     128 
     129      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition 
     130 
     131      ! gccd = matrix . gcdes 
     132      DO jj = 2, jpjm1 
     133         DO ji = fs_2, fs_jpim1   ! vector opt. 
     134            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   & 
     135               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   & 
     136               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   ) 
     137         END DO 
     138      END DO  
     139 
     140      ! alph = (gcr,gcr)/(gcdes,gccd) 
     141      radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  ) 
     142      IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain 
     143      alph = rnorme /radd 
     144 
     145      ! gcx = gcx + alph * gcdes 
     146      ! gcr = gcr - alph * gccd 
     147      DO jj = 2, jpjm1 
     148         DO ji = fs_2, fs_jpim1   ! vector opt. 
     149            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) ) 
     150            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) ) 
     151         END DO 
     152      END DO 
     153 
     154      ! Algorithm wtih Eijkhout rearrangement 
     155      ! ------------------------------------- 
     156         
    96157      !                                                !================ 
    97158      DO jn = 1, nmax                                  ! Iterative loop 
    98159         !                                             !================ 
    99160 
    100          IF( jn == 1 ) THEN           ! Initialization of the algorithm 
    101  
    102             CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition 
    103  
    104             ! gcr   = gcb-a.gcx 
    105             ! gcdes = gsr 
    106             DO jj = 2, jpjm1 
    107                DO ji = fs_2, fs_jpim1   ! vector opt. 
    108                   zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   & 
    109                      &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   & 
    110                      &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   & 
    111                      &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   & 
    112                      &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   ) 
    113                   gcr  (ji,jj) = zgcad 
    114                   gcdes(ji,jj) = zgcad 
    115                END DO 
    116             END DO 
    117              
    118             rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  ) 
    119             IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain 
    120             rr = rnorme 
    121  
    122          ENDIF 
    123          
    124          !                             ! Algorithm 
    125          
    126          CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition 
    127          
    128          ! ... gccd = matrix . gcdes 
     161         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition 
     162 
     163         ! zgcr = matrix . gcr 
    129164         DO jj = 2, jpjm1 
    130165            DO ji = fs_2, fs_jpim1   ! vector opt. 
    131                gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   & 
    132                   &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   & 
    133                   &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   ) 
     166               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   & 
     167                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   & 
     168                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   ) 
    134169            END DO 
    135170         END DO 
    136171  
    137          ! alph = (gcr,gcr)/(gcdes,gccd) 
    138          radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  ) 
    139          IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain 
    140          alph = rr / radd 
    141           
    142          ! gcx = gcx + alph * gcdes 
    143          ! gcr = gcr - alph * gccd 
    144          DO jj = 2, jpjm1 
    145             DO ji = fs_2, fs_jpim1   ! vector opt. 
    146                gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) ) 
    147                gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) ) 
    148             END DO 
    149          END DO 
    150          
    151172         ! rnorme = (gcr,gcr) 
    152          rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  ) 
    153          IF( lk_mpp )   CALL  mpp_sum( rnorme )   ! sum over the global domain 
    154          
     173         rr = rnorme 
     174         zsum(1) = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  ) 
     175 
     176         ! zgcad = (zgcr,gcr)  
     177         zsum(2) = SUM( gcr(2:jpim1,2:jpjm1) * gcdmat(2:jpim1,2:jpjm1) * zgcr(2:jpim1,2:jpjm1) ) 
     178 
     179         IF( lk_mpp )   CALL mpp_sum( zsum, 2 )   ! sum over the global domain 
     180         rnorme = zsum(1)   
     181         zgcad  = zsum(2) 
     182 
    155183         ! test of convergence 
    156184         IF( rnorme < epsr .OR. jn == nmax ) THEN 
     
    159187            ncut = 999 
    160188         ENDIF 
    161          
     189 
    162190         ! beta = (rk+1,rk+1)/(rk,rk) 
    163191         beta = rnorme / rr 
    164          rr   = rnorme 
    165  
     192         radd = zgcad - beta*beta*radd 
     193         alph = rnorme / radd 
     194 
     195         ! gcx = gcx + alph * gcdes 
     196         ! gcr = gcr - alph * gccd 
     197         DO jj = 2, jpjm1 
     198            DO ji = fs_2, fs_jpim1   ! vector opt. 
     199               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj)  
     200               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj)  
     201               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj)  
     202               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj)  
     203            END DO 
     204         END DO 
     205         
    166206         ! indicator of non-convergence or explosion 
    167207         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2 
    168208         IF( ncut == 999 ) GOTO 999 
    169209 
    170          ! gcdes = gcr + beta * gcdes 
    171          DO jj = 2, jpjm1 
    172             DO ji = fs_2, fs_jpim1   ! vector opt. 
    173                gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) ) 
    174             END DO 
    175          END DO 
    176          
    177210         !                                             !================ 
    178211      END DO                                           !    End Loop 
Note: See TracChangeset for help on using the changeset viewer.