Back to home page

MITgcm
 
 

    

File indexing completed on 2021-08-12 05:12:18 UTC

view on githubraw file Latest commit 0320e252 on 2021-08-11 16:08:52 UTC
66e5267b65 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002 
                0003 C--   File seaice_fgmres.F: seaice fgmres dynamical (linear) solver S/R:
                0004 C--   Contents
                0005 C--   o SEAICE_MAP2VEC
ae95ff41d3 Jean*0006 C--   o SEAICE_MAP_RS2VEC
66e5267b65 Mart*0007 C--   o SEAICE_FGMRES
ae37ac4a14 Mart*0008 C--   o SEAICE_SCALPROD
66e5267b65 Mart*0009 
                0010 CBOP
                0011 C     !ROUTINE: SEAICE_MAP2VEC
                0012 C     !INTERFACE:
                0013 
                0014       SUBROUTINE SEAICE_MAP2VEC(
372e6d06fe Jean*0015      I     n,
                0016      O     xfld2d, yfld2d,
                0017      U     vector,
66e5267b65 Mart*0018      I     map2vec, myThid )
                0019 
                0020 C     !DESCRIPTION: \bv
                0021 C     *==========================================================*
                0022 C     | SUBROUTINE SEAICE_MAP2VEC
                0023 C     | o maps 2 2D-fields to vector and back
                0024 C     *==========================================================*
                0025 C     | written by Martin Losch, Oct 2012
                0026 C     *==========================================================*
                0027 C     \ev
                0028 
                0029 C     !USES:
                0030       IMPLICIT NONE
                0031 
                0032 C     === Global variables ===
                0033 #include "SIZE.h"
                0034 #include "EEPARAMS.h"
                0035 C     === Routine arguments ===
                0036       INTEGER n
                0037       LOGICAL map2vec
                0038       INTEGER myThid
                0039       _RL xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0040       _RL yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
dd78fdcb7e Mart*0041       _RL vector (n,nSx,nSy)
7871461369 Mart*0042 #if (defined SEAICE_ALLOW_JFNK) || (defined SEAICE_ALLOW_KRYLOV)
66e5267b65 Mart*0043 C     === local variables ===
                0044       INTEGER I, J, bi, bj
dd78fdcb7e Mart*0045       INTEGER ii, jj, m
66e5267b65 Mart*0046 CEOP
372e6d06fe Jean*0047 
66e5267b65 Mart*0048       m = n/2
dd78fdcb7e Mart*0049       DO bj=myByLo(myThid),myByHi(myThid)
                0050        DO bi=myBxLo(myThid),myBxHi(myThid)
93a55af7df Mart*0051 #ifdef SEAICE_JFNK_MAP_REORDER
                0052         ii = 0
                0053         IF ( map2vec ) THEN
                0054          DO J=1,sNy
                0055           jj = 2*sNx*(J-1)
                0056           DO I=1,sNx
                0057            ii = jj + 2*I
                0058            vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
                0059            vector(ii,  bi,bj) = yfld2d(I,J,bi,bj)
                0060           ENDDO
                0061          ENDDO
                0062         ELSE
                0063          DO J=1,sNy
                0064           jj = 2*sNx*(J-1)
                0065           DO I=1,sNx
                0066            ii = jj + 2*I
                0067            xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
                0068            yfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
                0069           ENDDO
                0070          ENDDO
                0071         ENDIF
                0072 #else
dd78fdcb7e Mart*0073         IF ( map2vec ) THEN
66e5267b65 Mart*0074          DO J=1,sNy
dd78fdcb7e Mart*0075           jj = sNx*(J-1)
66e5267b65 Mart*0076           DO I=1,sNx
                0077            ii = jj + I
dd78fdcb7e Mart*0078            vector(ii,  bi,bj) = xfld2d(I,J,bi,bj)
                0079            vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
66e5267b65 Mart*0080           ENDDO
                0081          ENDDO
dd78fdcb7e Mart*0082         ELSE
66e5267b65 Mart*0083          DO J=1,sNy
dd78fdcb7e Mart*0084           jj = sNx*(J-1)
66e5267b65 Mart*0085           DO I=1,sNx
                0086            ii = jj + I
dd78fdcb7e Mart*0087            xfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
                0088            yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
66e5267b65 Mart*0089           ENDDO
                0090          ENDDO
dd78fdcb7e Mart*0091         ENDIF
93a55af7df Mart*0092 #endif /* SEAICE_JFNK_MAP_REORDER */
dd78fdcb7e Mart*0093 C     bi,bj-loops
372e6d06fe Jean*0094        ENDDO
dd78fdcb7e Mart*0095       ENDDO
66e5267b65 Mart*0096 
7871461369 Mart*0097 #endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
66e5267b65 Mart*0098       RETURN
                0099       END
                0100 
                0101 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0102 CBOP
ae95ff41d3 Jean*0103 C     !ROUTINE: SEAICE_MAP_RS2VEC
                0104 C     !INTERFACE:
                0105 
                0106       SUBROUTINE SEAICE_MAP_RS2VEC(
                0107      I     n,
                0108      O     xfld2d, yfld2d,
                0109      U     vector,
                0110      I     map2vec, myThid )
                0111 
                0112 C     !DESCRIPTION: \bv
                0113 C     *==========================================================*
                0114 C     | SUBROUTINE SEAICE_MAP_RS2VEC
                0115 C     | o maps 2 2D-RS-fields to vector and back
                0116 C     *==========================================================*
                0117 C     | written by Martin Losch, Oct 2012
                0118 C     *==========================================================*
                0119 C     \ev
                0120 
                0121 C     !USES:
                0122       IMPLICIT NONE
                0123 
                0124 C     === Global variables ===
                0125 #include "SIZE.h"
                0126 #include "EEPARAMS.h"
                0127 C     === Routine arguments ===
                0128       INTEGER n
                0129       LOGICAL map2vec
                0130       INTEGER myThid
                0131       _RS xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0132       _RS yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0133       _RL vector (n,nSx,nSy)
7871461369 Mart*0134 #if (defined SEAICE_ALLOW_JFNK) || (defined SEAICE_ALLOW_KRYLOV)
ae95ff41d3 Jean*0135 C     === local variables ===
                0136       INTEGER I, J, bi, bj
                0137       INTEGER ii, jj, m
                0138 CEOP
                0139 
                0140       m = n/2
                0141       DO bj=myByLo(myThid),myByHi(myThid)
                0142        DO bi=myBxLo(myThid),myBxHi(myThid)
93a55af7df Mart*0143 #ifdef SEAICE_JFNK_MAP_REORDER
                0144         ii = 0
                0145         IF ( map2vec ) THEN
                0146          DO J=1,sNy
                0147           jj = 2*sNx*(J-1)
                0148           DO I=1,sNx
                0149            ii = jj + 2*I
                0150            vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
                0151            vector(ii,  bi,bj) = yfld2d(I,J,bi,bj)
                0152           ENDDO
                0153          ENDDO
                0154         ELSE
                0155          DO J=1,sNy
                0156           jj = 2*sNx*(J-1)
                0157           DO I=1,sNx
                0158            ii = jj + 2*I
                0159            xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
                0160            yfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
                0161           ENDDO
                0162          ENDDO
                0163         ENDIF
                0164 #else
ae95ff41d3 Jean*0165         IF ( map2vec ) THEN
                0166          DO J=1,sNy
                0167           jj = sNx*(J-1)
                0168           DO I=1,sNx
                0169            ii = jj + I
                0170            vector(ii,  bi,bj) = xfld2d(I,J,bi,bj)
                0171            vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
                0172           ENDDO
                0173          ENDDO
                0174         ELSE
                0175          DO J=1,sNy
                0176           jj = sNx*(J-1)
                0177           DO I=1,sNx
                0178            ii = jj + I
                0179            xfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
                0180            yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
                0181           ENDDO
                0182          ENDDO
                0183         ENDIF
93a55af7df Mart*0184 #endif /* SEAICE_JFNK_MAP_REORDER */
ae95ff41d3 Jean*0185 C     bi,bj-loops
                0186        ENDDO
                0187       ENDDO
                0188 
7871461369 Mart*0189 #endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
ae95ff41d3 Jean*0190       RETURN
                0191       END
                0192 
                0193 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0194 CBOP
66e5267b65 Mart*0195 C     !ROUTINE: SEAICE_FGMRES
                0196 C     !INTERFACE:
ae37ac4a14 Mart*0197       SUBROUTINE SEAICE_FGMRES (
                0198      I     n,im,rhs,
                0199      U     sol,i,its,vv,w,wk1,wk2,
                0200      I     eps,maxits,iout,
                0201      U     icode,
                0202      I     myThid )
66e5267b65 Mart*0203 
                0204 C-----------------------------------------------------------------------
                0205 C mlosch Oct 2012: modified the routine further to be compliant with
                0206 C MITgcm standards:
                0207 C f90 -> F
                0208 C !-comment -> C-comment
ae37ac4a14 Mart*0209 C add its to list of arguments
66e5267b65 Mart*0210 C double precision -> _RL
                0211 C implicit none
372e6d06fe Jean*0212 C
66e5267b65 Mart*0213 C jfl Dec 1st 2006. We modified the routine so that it is double precison.
                0214 C Here are the modifications:
372e6d06fe Jean*0215 C 1) implicit real (a-h,o-z) becomes implicit real*8 (a-h,o-z)
66e5267b65 Mart*0216 C 2) real bocomes real*8
                0217 C 3) subroutine scopy.f has been changed for dcopy.f
                0218 C 4) subroutine saxpy.f has been changed for daxpy.f
                0219 C 5) function sdot.f has been changed for ddot.f
                0220 C 6) 1e-08 becomes 1d-08
                0221 C
372e6d06fe Jean*0222 C Be careful with the dcopy, daxpy and ddot code...there is a slight
66e5267b65 Mart*0223 C difference with the single precision versions (scopy, saxpy and sdot).
                0224 C In the single precision versions, the array are declared sightly differently.
                0225 C It is written for single precision:
                0226 C
                0227 C modified 12/3/93, array(1) declarations changed to array(*)
                0228 C-----------------------------------------------------------------------
                0229 
                0230       implicit none
ae37ac4a14 Mart*0231 C     === Global variables ===
dd78fdcb7e Mart*0232 #include "SIZE.h"
                0233 #include "EEPARAMS.h"
66e5267b65 Mart*0234       integer myThid
7871461369 Mart*0235       integer i, n, im, its, maxits, iout, icode
dd78fdcb7e Mart*0236       _RL rhs(n,nSx,nSy), sol(n,nSx,nSy)
                0237       _RL vv(n,im+1,nSx,nSy), w(n,im,nSx,nSy)
                0238       _RL wk1(n,nSx,nSy), wk2(n,nSx,nSy), eps
66e5267b65 Mart*0239 C-----------------------------------------------------------------------
372e6d06fe Jean*0240 C flexible GMRES routine. This is a version of GMRES which allows a
                0241 C a variable preconditioner. Implemented with a reverse communication
66e5267b65 Mart*0242 C protocole for flexibility -
372e6d06fe Jean*0243 C DISTRIBUTED VERSION (USES DISTDOT FOR DDOT)
                0244 C explicit (exact) residual norms for restarts
66e5267b65 Mart*0245 C written by Y. Saad, modified by A. Malevsky, version February 1, 1995
                0246 C-----------------------------------------------------------------------
372e6d06fe Jean*0247 C This Is A Reverse Communication Implementation.
                0248 C-------------------------------------------------
66e5267b65 Mart*0249 C USAGE: (see also comments for icode below). FGMRES
                0250 C should be put in a loop and the loop should be active for as
                0251 C long as icode is not equal to 0. On return fgmres will
                0252 C    1) either be requesting the new preconditioned vector applied
372e6d06fe Jean*0253 C       to wk1 in case icode.eq.1 (result should be put in wk2)
66e5267b65 Mart*0254 C    2) or be requesting the product of A applied to the vector wk1
372e6d06fe Jean*0255 C       in case icode.eq.2 (result should be put in wk2)
                0256 C    3) or be terminated in case icode .eq. 0.
66e5267b65 Mart*0257 C on entry always set icode = 0. So icode should be set back to zero
                0258 C upon convergence.
                0259 C-----------------------------------------------------------------------
372e6d06fe Jean*0260 C Here is a typical way of running fgmres:
66e5267b65 Mart*0261 C
                0262 C      icode = 0
                0263 C 1    continue
ae37ac4a14 Mart*0264 C      call fgmres (n,im,rhs,sol,i,vv,w,wk1,wk2,eps,maxits,iout,
                0265 C     &             icode,its,mythid)
66e5267b65 Mart*0266 C
                0267 C      if (icode .eq. 1) then
c04a128844 Patr*0268 C         call  precon(n, wk1, wk2)    <--- user variable preconditioning
66e5267b65 Mart*0269 C         goto 1
                0270 C      else if (icode .ge. 2) then
372e6d06fe Jean*0271 C         call  matvec (n,wk1, wk2)    <--- user matrix vector product.
66e5267b65 Mart*0272 C         goto 1
372e6d06fe Jean*0273 C      else
                0274 C         ----- done ----
66e5267b65 Mart*0275 C         .........
                0276 C-----------------------------------------------------------------------
372e6d06fe Jean*0277 C list of parameters
                0278 C-------------------
66e5267b65 Mart*0279 C
                0280 C n     == integer. the dimension of the problem
                0281 C im    == size of Krylov subspace:  should not exceed 50 in this
                0282 C          version (can be reset in code. looking at comment below)
                0283 C rhs   == vector of length n containing the right hand side
                0284 C sol   == initial guess on input, approximate solution on output
                0285 C vv    == work space of size n x (im+1)
372e6d06fe Jean*0286 C w     == work space of length n x im
66e5267b65 Mart*0287 C wk1,
                0288 C wk2,  == two work vectors of length n each used for the reverse
                0289 C          communication protocole. When on return (icode .ne. 1)
                0290 C          the user should call fgmres again with wk2 = precon * wk1
                0291 C          and icode untouched. When icode.eq.1 then it means that
                0292 C          convergence has taken place.
372e6d06fe Jean*0293 C
66e5267b65 Mart*0294 C eps   == tolerance for stopping criterion. process is stopped
                0295 C          as soon as ( ||.|| is the euclidean norm):
                0296 C          || current residual||/||initial residual|| <= eps
                0297 C
                0298 C maxits== maximum number of iterations allowed
                0299 C
ae37ac4a14 Mart*0300 C i     == internal iteration counter, updated in this routine
                0301 C its   == current (Krylov) iteration counter, updated in this routine
                0302 C
66e5267b65 Mart*0303 C iout  == output unit number number for printing intermediate results
                0304 C          if (iout .le. 0) no statistics are printed.
372e6d06fe Jean*0305 C
66e5267b65 Mart*0306 C icode = integer. indicator for the reverse communication protocole.
                0307 C         ON ENTRY : icode should be set to icode = 0.
372e6d06fe Jean*0308 C         ON RETURN:
66e5267b65 Mart*0309 C       * icode .eq. 1 value means that fgmres has not finished
                0310 C         and that it is requesting a preconditioned vector before
                0311 C         continuing. The user must compute M**(-1) wk1, where M is
                0312 C         the preconditioing  matrix (may vary at each call) and wk1 is
372e6d06fe Jean*0313 C         the vector as provided by fgmres upun return, and put the
66e5267b65 Mart*0314 C         result in wk2. Then fgmres must be called again without
372e6d06fe Jean*0315 C         changing any other argument.
66e5267b65 Mart*0316 C       * icode .eq. 2 value means that fgmres has not finished
                0317 C         and that it is requesting a matrix vector product before
                0318 C         continuing. The user must compute  A * wk1, where A is the
372e6d06fe Jean*0319 C         coefficient  matrix and wk1 is the vector provided by
66e5267b65 Mart*0320 C         upon return. The result of the operation is to be put in
                0321 C         the vector wk2. Then fgmres must be called again without
372e6d06fe Jean*0322 C         changing any other argument.
                0323 C       * icode .eq. 0 means that fgmres has finished and sol contains
66e5267b65 Mart*0324 C         the approximate solution.
                0325 C         comment: typically fgmres must be implemented in a loop
372e6d06fe Jean*0326 C         with fgmres being called as long icode is returned with
                0327 C         a value .ne. 0.
66e5267b65 Mart*0328 C-----------------------------------------------------------------------
7871461369 Mart*0329 #if (defined SEAICE_ALLOW_JFNK) || (defined SEAICE_ALLOW_KRYLOV)
66e5267b65 Mart*0330 C     local variables -- !jfl modif
372e6d06fe Jean*0331       integer imax
66e5267b65 Mart*0332       parameter ( imax = 50 )
b3609d72bf Mart*0333       _RL hh(4*imax+1,4*imax,MAX_NO_THREADS)
                0334       _RL c(4*imax,MAX_NO_THREADS),s(4*imax,MAX_NO_THREADS)
                0335       _RL rs(4*imax+1,MAX_NO_THREADS),t,ro
66e5267b65 Mart*0336 C-------------------------------------------------------------
                0337 C     arnoldi size should not exceed 50 in this version..
                0338 C-------------------------------------------------------------
b3609d72bf Mart*0339       integer i1Thid(MAX_NO_THREADS)
7871461369 Mart*0340       integer i1, ii, j, jj, k, k1!, n1
dd78fdcb7e Mart*0341       integer bi, bj
2084da270c Mart*0342       _RL r0, gam, eps1(MAX_NO_THREADS)
dd78fdcb7e Mart*0343       CHARACTER*(MAX_LEN_MBUF) msgBuf
66e5267b65 Mart*0344 
                0345 CEOP
b3609d72bf Mart*0346 C     local common blocks to replace a save statement in the original code
                0347       COMMON /SEAICE_FMRES_LOC_I/ i1Thid
80dd7a47d1 Mart*0348       COMMON /SEAICE_FMRES_LOC_RL/ hh, c, s, rs, eps1
                0349       _RL epsmac
                0350       parameter ( epsmac =  1.d-16 )
372e6d06fe Jean*0351 C
                0352 C     computed goto
                0353 C
66e5267b65 Mart*0354       if ( im .gt. imax ) stop 'size of krylov space > 50'
                0355       goto (100,200,300,11) icode +1
                0356  100  continue
                0357       its = 0
                0358 C-------------------------------------------------------------
                0359 C     **  outer loop starts here..
                0360 C--------------compute initial residual vector --------------
                0361 C 10   continue
dd78fdcb7e Mart*0362       do bj=myByLo(myThid),myByHi(myThid)
                0363        do bi=myBxLo(myThid),myBxHi(myThid)
93a55af7df Mart*0364         do j=1,n
                0365          wk1(j,bi,bj)=sol(j,bi,bj)
dd78fdcb7e Mart*0366         enddo
                0367        enddo
66e5267b65 Mart*0368       enddo
                0369       icode = 3
372e6d06fe Jean*0370       RETURN
66e5267b65 Mart*0371  11   continue
dd78fdcb7e Mart*0372       do bj=myByLo(myThid),myByHi(myThid)
                0373        do bi=myBxLo(myThid),myBxHi(myThid)
                0374         do j=1,n
                0375          vv(j,1,bi,bj) = rhs(j,bi,bj) - wk2(j,bi,bj)
                0376         enddo
                0377        enddo
66e5267b65 Mart*0378       enddo
2899eda9c1 Mart*0379  20   continue
                0380       call SEAICE_SCALPROD(n, im+1, 1, 1, vv, vv, ro, myThid)
66e5267b65 Mart*0381       ro = sqrt(ro)
438b73e6d6 Mart*0382       if (ro .eq. 0.0 _d 0) goto 999
                0383       t = 1.0 _d 0/ ro
dd78fdcb7e Mart*0384       do bj=myByLo(myThid),myByHi(myThid)
                0385        do bi=myBxLo(myThid),myBxHi(myThid)
                0386         do j=1, n
                0387          vv(j,1,bi,bj) = vv(j,1,bi,bj)*t
                0388         enddo
                0389        enddo
66e5267b65 Mart*0390       enddo
2084da270c Mart*0391       if (its .eq. 0) eps1(myThid)=eps
dd78fdcb7e Mart*0392 C     not sure what this is, r0 is never used again
66e5267b65 Mart*0393       if (its .eq. 0) r0 = ro
dd78fdcb7e Mart*0394       if (iout .gt. 0) then
ab70c5fd89 Jean*0395        _BEGIN_MASTER( myThid )
dd78fdcb7e Mart*0396        write(msgBuf, 199) its, ro
                0397        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0398      &      SQUEEZE_RIGHT, myThid )
66e5267b65 Mart*0399 C           print *,'chau',its, ro !write(iout, 199) its, ro
dd78fdcb7e Mart*0400        _END_MASTER( myThid )
                0401       endif
372e6d06fe Jean*0402 C
66e5267b65 Mart*0403 C     initialize 1-st term  of rhs of hessenberg system..
372e6d06fe Jean*0404 C
b3609d72bf Mart*0405       rs(1,myThid) = ro
66e5267b65 Mart*0406       i = 0
2899eda9c1 Mart*0407  4    continue
                0408       i=i+1
66e5267b65 Mart*0409       its = its + 1
b3609d72bf Mart*0410       i1Thid(myThid) = i + 1
dd78fdcb7e Mart*0411       do bj=myByLo(myThid),myByHi(myThid)
                0412        do bi=myBxLo(myThid),myBxHi(myThid)
                0413         do k=1, n
                0414          wk1(k,bi,bj) = vv(k,i,bi,bj)
                0415         enddo
                0416        enddo
66e5267b65 Mart*0417       enddo
372e6d06fe Jean*0418 C
66e5267b65 Mart*0419 C     return
372e6d06fe Jean*0420 C
66e5267b65 Mart*0421       icode = 1
372e6d06fe Jean*0422       RETURN
66e5267b65 Mart*0423  200  continue
dd78fdcb7e Mart*0424       do bj=myByLo(myThid),myByHi(myThid)
                0425        do bi=myBxLo(myThid),myBxHi(myThid)
                0426         do k=1, n
                0427          w(k,i,bi,bj) = wk2(k,bi,bj)
                0428         enddo
                0429        enddo
66e5267b65 Mart*0430       enddo
372e6d06fe Jean*0431 C
66e5267b65 Mart*0432 C     call matvec operation
372e6d06fe Jean*0433 C
dd78fdcb7e Mart*0434       do bj=myByLo(myThid),myByHi(myThid)
                0435        do bi=myBxLo(myThid),myBxHi(myThid)
                0436         do k=1,n
                0437          wk1(k,bi,bj)=wk2(k,bi,bj)
                0438         enddo
                0439        enddo
66e5267b65 Mart*0440       enddo
                0441 C
                0442 C     return
372e6d06fe Jean*0443 C
2899eda9c1 Mart*0444       icode = 2
372e6d06fe Jean*0445       RETURN
66e5267b65 Mart*0446  300  continue
372e6d06fe Jean*0447 C
66e5267b65 Mart*0448 C     first call to ope corresponds to intialization goto back to 11.
372e6d06fe Jean*0449 C
66e5267b65 Mart*0450 C      if (icode .eq. 3) goto 11
b3609d72bf Mart*0451       i1 = i1Thid(myThid)
dd78fdcb7e Mart*0452       do bj=myByLo(myThid),myByHi(myThid)
                0453        do bi=myBxLo(myThid),myBxHi(myThid)
                0454         do k=1,n
                0455          vv(k,i1,bi,bj)=wk2(k,bi,bj)
                0456         enddo
                0457        enddo
66e5267b65 Mart*0458       enddo
372e6d06fe Jean*0459 C
66e5267b65 Mart*0460 C     modified gram - schmidt...
372e6d06fe Jean*0461 C
66e5267b65 Mart*0462       do j=1, i
ae37ac4a14 Mart*0463        call SEAICE_SCALPROD(n, im+1, j, i1, vv, vv, t, myThid)
b3609d72bf Mart*0464        hh(j,i,myThid) = t
80dd7a47d1 Mart*0465        do bj=myByLo(myThid),myByHi(myThid)
                0466         do bi=myBxLo(myThid),myBxHi(myThid)
66e5267b65 Mart*0467          do k=1,n
dd78fdcb7e Mart*0468           vv(k,i1,bi,bj) = vv(k,i1,bi,bj) - t*vv(k,j,bi,bj)
66e5267b65 Mart*0469          enddo
dd78fdcb7e Mart*0470         enddo
                0471        enddo
66e5267b65 Mart*0472       enddo
ae37ac4a14 Mart*0473       call SEAICE_SCALPROD(n, im+1, i1, i1, vv, vv, t, myThid)
66e5267b65 Mart*0474       t = sqrt(t)
b3609d72bf Mart*0475       hh(i1,i,myThid) = t
438b73e6d6 Mart*0476       if (t .ne. 0.0 _d 0) then
                0477        t = 1.0 _d 0 / t
dd78fdcb7e Mart*0478        do bj=myByLo(myThid),myByHi(myThid)
                0479         do bi=myBxLo(myThid),myBxHi(myThid)
                0480          do k=1,n
                0481           vv(k,i1,bi,bj) = vv(k,i1,bi,bj)*t
                0482          enddo
                0483         enddo
                0484        enddo
                0485       endif
372e6d06fe Jean*0486 C
80dd7a47d1 Mart*0487 C     done with modified gram schmidt and arnoldi step.
66e5267b65 Mart*0488 C     now  update factorization of hh
372e6d06fe Jean*0489 C
b3609d72bf Mart*0490       if (i .ne. 1) then
372e6d06fe Jean*0491 C
66e5267b65 Mart*0492 C     perfrom previous transformations  on i-th column of h
372e6d06fe Jean*0493 C
b3609d72bf Mart*0494        do k=2,i
                0495         k1 = k-1
                0496         t = hh(k1,i,myThid)
                0497         hh(k1,i,myThid) = c(k1,myThid)*t+s(k1,myThid)*hh(k,i,myThid)
                0498         hh(k,i,myThid) = -s(k1,myThid)*t+c(k1,myThid)*hh(k,i,myThid)
                0499        enddo
                0500       endif
                0501       gam = sqrt(hh(i,i,myThid)**2 + hh(i1,i,myThid)**2)
                0502       if (gam .eq. 0.0 _d 0) gam = epsmac
66e5267b65 Mart*0503 C-----------#determine next plane rotation  #-------------------
b3609d72bf Mart*0504       c(i,myThid)   = hh(i, i,myThid)/gam
                0505       s(i,myThid)   = hh(i1,i,myThid)/gam
438b73e6d6 Mart*0506 C     numerically more stable Givens rotation, but the results
                0507 C     are not better
b3609d72bf Mart*0508 CML      c(i,myThid)=1. _d 0
                0509 CML      s(i,myThid)=0. _d 0
                0510 CML      if ( abs(hh(i1(myThid,myThid),i)) .gt. 0.0 _d 0) then
                0511 CML       if ( abs(hh(i1,i,myThid)) .gt. abs(hh(i,i,myThid)) ) then
                0512 CML        gam = hh(i,i,myThid)/hh(i1,i,myThid)
                0513 CML        s(i,myThid) = 1./sqrt(1.+gam*gam)
                0514 CML        c(i,myThid) = s(i,myThid)*gam
                0515 CML       else
                0516 CML        gam = hh(i1,i,myThid)/hh(i,i,myThid)
                0517 CML        c(i,myThid) = 1./sqrt(1.+gam*gam)
                0518 CML        s(i,myThid) = c(i,myThid)*gam
                0519 CML       endif
                0520 CML      endif
                0521       rs(i1,myThid) = -s(i,myThid)*rs(i,myThid)
                0522       rs(i ,myThid) =  c(i,myThid)*rs(i,myThid)
372e6d06fe Jean*0523 C
2899eda9c1 Mart*0524 C     determine res. norm. and test for convergence
372e6d06fe Jean*0525 C
0320e25227 Mart*0526       hh(i,i,myThid) = c(i,myThid)*hh(i, i,myThid)
b3609d72bf Mart*0527      &     +           s(i,myThid)*hh(i1,i,myThid)
                0528       ro = abs(rs(i1,myThid))
dd78fdcb7e Mart*0529       if (iout .gt. 0) then
ab70c5fd89 Jean*0530        _BEGIN_MASTER( myThid )
dd78fdcb7e Mart*0531        write(msgBuf, 199) its, ro
                0532        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0533      &      SQUEEZE_RIGHT, myThid )
                0534        _END_MASTER( myThid )
                0535       endif
a8d81a425a Mart*0536 CML   The original code uses only im as the maximum allowed iterations, but
                0537 CML   I want no more than maxits iterations, so we change the code here.
2084da270c Mart*0538 CML      if (i .lt. im .and. (ro .gt. eps1(myThid)))  goto 4
0320e25227 Mart*0539       if (i .lt. im .and. its .lt. maxits
2084da270c Mart*0540      &     .and. (ro .gt. eps1(myThid))) goto 4
372e6d06fe Jean*0541 C
66e5267b65 Mart*0542 C     now compute solution. first solve upper triangular system.
372e6d06fe Jean*0543 C
b3609d72bf Mart*0544       rs(i,myThid) = rs(i,myThid)/hh(i,i,myThid)
                0545       do ii=2,i
                0546        k=i-ii+1
                0547        k1 = k+1
                0548        t=rs(k,myThid)
                0549        do j=k1,i
                0550         t = t-hh(k,j,myThid)*rs(j,myThid)
2899eda9c1 Mart*0551        enddo
b3609d72bf Mart*0552        rs(k,myThid) = t/hh(k,k,myThid)
66e5267b65 Mart*0553       enddo
372e6d06fe Jean*0554 C
66e5267b65 Mart*0555 C     done with back substitution..
                0556 C     now form linear combination to get solution
372e6d06fe Jean*0557 C
66e5267b65 Mart*0558       do j=1, i
b3609d72bf Mart*0559        t = rs(j,myThid)
dd78fdcb7e Mart*0560        do bj=myByLo(myThid),myByHi(myThid)
                0561         do bi=myBxLo(myThid),myBxHi(myThid)
                0562          do k=1,n
b3609d72bf Mart*0563           sol(k,bi,bj) = sol(k,bi,bj) + t*w(k,j,bi,bj)
dd78fdcb7e Mart*0564          enddo
                0565         enddo
66e5267b65 Mart*0566        enddo
                0567       enddo
372e6d06fe Jean*0568 C
                0569 C     test for return
                0570 C
2084da270c Mart*0571       if (ro .le. eps1(myThid) .or. its .ge. maxits) goto 999
372e6d06fe Jean*0572 C
66e5267b65 Mart*0573 C     else compute residual vector and continue..
372e6d06fe Jean*0574 C
66e5267b65 Mart*0575 C       goto 10
                0576 
b3609d72bf Mart*0577       i1 = i1Thid(myThid)
                0578       do j=1,i
                0579        jj = i1-j+1
                0580        rs(jj-1,myThid) = -s(jj-1,myThid)*rs(jj,myThid)
                0581        rs(jj  ,myThid) =  c(jj-1,myThid)*rs(jj,myThid)
                0582       enddo
                0583       do j=1,i1
                0584        t = rs(j,myThid)
                0585        if (j .eq. 1)  t = t - 1.0 _d 0
                0586        do bj=myByLo(myThid),myByHi(myThid)
                0587         do bi=myBxLo(myThid),myBxHi(myThid)
2899eda9c1 Mart*0588          do k=1,n
                0589           vv(k,1,bi,bj) = vv(k,1,bi,bj) + t*vv(k,j,bi,bj)
dd78fdcb7e Mart*0590          enddo
66e5267b65 Mart*0591         enddo
2899eda9c1 Mart*0592        enddo
66e5267b65 Mart*0593       enddo
372e6d06fe Jean*0594 C
66e5267b65 Mart*0595 C     restart outer loop.
372e6d06fe Jean*0596 C
66e5267b65 Mart*0597       goto 20
                0598  999  icode = 0
                0599 
ae37ac4a14 Mart*0600  199  format(' SEAICE_FGMRES: its =', i4, ' res. norm =', d26.16)
372e6d06fe Jean*0601 C
7871461369 Mart*0602 #endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
372e6d06fe Jean*0603       RETURN
                0604 C-----end-of-fgmres-----------------------------------------------------
66e5267b65 Mart*0605 C-----------------------------------------------------------------------
372e6d06fe Jean*0606       END
66e5267b65 Mart*0607 
                0608 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0609 CBOP
ae37ac4a14 Mart*0610 C     !ROUTINE: SEAICE_SCALPROD
66e5267b65 Mart*0611 C     !INTERFACE:
                0612 
ae37ac4a14 Mart*0613       subroutine SEAICE_SCALPROD(n,im,i1,i2,dx,dy,t,myThid)
66e5267b65 Mart*0614 
                0615 C     forms the dot product of two vectors.
                0616 C     uses unrolled loops for increments equal to one.
                0617 C     jack dongarra, linpack, 3/11/78.
dd78fdcb7e Mart*0618 C     ML: code stolen from BLAS-ddot and adapted for parallel applications
372e6d06fe Jean*0619 
66e5267b65 Mart*0620       implicit none
                0621 #include "SIZE.h"
                0622 #include "EEPARAMS.h"
                0623 #include "EESUPPORT.h"
438b73e6d6 Mart*0624 #include "SEAICE_SIZE.h"
                0625 #include "SEAICE.h"
dd78fdcb7e Mart*0626       integer n, im, i1, i2
ae37ac4a14 Mart*0627       _RL dx(n,im,nSx,nSy),dy(n,im,nSx,nSy)
dd78fdcb7e Mart*0628       _RL t
372e6d06fe Jean*0629       integer myThid
7871461369 Mart*0630 #if (defined SEAICE_ALLOW_JFNK) || (defined SEAICE_ALLOW_KRYLOV)
dd78fdcb7e Mart*0631 C     local arrays
                0632       _RL dtemp(nSx,nSy)
                0633       integer i,m,mp1,bi,bj
66e5267b65 Mart*0634 CEOP
ab70c5fd89 Jean*0635 
66e5267b65 Mart*0636       m = mod(n,5)
ab70c5fd89 Jean*0637       mp1 = m + 1
66e5267b65 Mart*0638       t     = 0. _d 0
ab70c5fd89 Jean*0639 c     if( m .eq. 0 ) go to 40
dd78fdcb7e Mart*0640       do bj=myByLo(myThid),myByHi(myThid)
                0641        do bi=myBxLo(myThid),myBxHi(myThid)
                0642         dtemp(bi,bj) = 0. _d 0
ab70c5fd89 Jean*0643         if ( m .ne. 0 ) then
                0644          do i = 1,m
                0645           dtemp(bi,bj) = dtemp(bi,bj) + dx(i,i1,bi,bj)*dy(i,i2,bi,bj)
e5c7689f4a Mart*0646      &         * scalarProductMetric(i,1,bi,bj)
ab70c5fd89 Jean*0647          enddo
                0648         endif
                0649         if ( n .ge. 5 ) then
                0650 c     if( n .lt. 5 ) go to 60
                0651 c40   mp1 = m + 1
                0652          do i = mp1,n,5
438b73e6d6 Mart*0653           dtemp(bi,bj) = dtemp(bi,bj)               +
                0654      &        dx(i,    i1,bi,bj)*dy(i,    i2,bi,bj)
                0655      &        * scalarProductMetric(i,    1, bi,bj) +
                0656      &        dx(i + 1,i1,bi,bj)*dy(i + 1,i2,bi,bj)
                0657      &        * scalarProductMetric(i + 1,1, bi,bj) +
                0658      &        dx(i + 2,i1,bi,bj)*dy(i + 2,i2,bi,bj)
                0659      &        * scalarProductMetric(i + 2,1, bi,bj) +
                0660      &        dx(i + 3,i1,bi,bj)*dy(i + 3,i2,bi,bj)
                0661      &        * scalarProductMetric(i + 3,1, bi,bj) +
dd78fdcb7e Mart*0662      &        dx(i + 4,i1,bi,bj)*dy(i + 4,i2,bi,bj)
438b73e6d6 Mart*0663      &        * scalarProductMetric(i + 4,1, bi,bj)
ab70c5fd89 Jean*0664          enddo
                0665 c60   continue
                0666         endif
dd78fdcb7e Mart*0667        enddo
66e5267b65 Mart*0668       enddo
dd78fdcb7e Mart*0669       CALL GLOBAL_SUM_TILE_RL( dtemp,t,myThid )
372e6d06fe Jean*0670 
7871461369 Mart*0671 #endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
66e5267b65 Mart*0672       RETURN
                0673       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!