C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_fgmres.F,v 1.22 2016/09/27 05:34:46 mlosch Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
C-- File seaice_fgmres.F: seaice fgmres dynamical (linear) solver S/R:
C-- Contents
C-- o SEAICE_MAP2VEC
C-- o SEAICE_MAP_RS2VEC
C-- o SEAICE_FGMRES
C-- o SEAICE_SCALPROD
CBOP
C !ROUTINE: SEAICE_MAP2VEC
C !INTERFACE:
SUBROUTINE SEAICE_MAP2VEC(
I n,
O xfld2d, yfld2d,
U vector,
I map2vec, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SEAICE_MAP2VEC
C | o maps 2 2D-fields to vector and back
C *==========================================================*
C | written by Martin Losch, Oct 2012
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
C === Routine arguments ===
INTEGER n
LOGICAL map2vec
INTEGER myThid
_RL xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vector (n,nSx,nSy)
#if (defined SEAICE_ALLOW_JFNK) (defined SEAICE_ALLOW_KRYLOV)
C === local variables ===
INTEGER I, J, bi, bj
INTEGER ii, jj, m
CEOP
m = n/2
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef SEAICE_JFNK_MAP_REORDER
ii = 0
IF ( map2vec ) THEN
DO J=1,sNy
jj = 2*sNx*(J-1)
DO I=1,sNx
ii = jj + 2*I
vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
vector(ii, bi,bj) = yfld2d(I,J,bi,bj)
ENDDO
ENDDO
ELSE
DO J=1,sNy
jj = 2*sNx*(J-1)
DO I=1,sNx
ii = jj + 2*I
xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
yfld2d(I,J,bi,bj) = vector(ii, bi,bj)
ENDDO
ENDDO
ENDIF
#else
IF ( map2vec ) THEN
DO J=1,sNy
jj = sNx*(J-1)
DO I=1,sNx
ii = jj + I
vector(ii, bi,bj) = xfld2d(I,J,bi,bj)
vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
ENDDO
ENDDO
ELSE
DO J=1,sNy
jj = sNx*(J-1)
DO I=1,sNx
ii = jj + I
xfld2d(I,J,bi,bj) = vector(ii, bi,bj)
yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* SEAICE_JFNK_MAP_REORDER */
C bi,bj-loops
ENDDO
ENDDO
#endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: SEAICE_MAP_RS2VEC
C !INTERFACE:
SUBROUTINE SEAICE_MAP_RS2VEC(
I n,
O xfld2d, yfld2d,
U vector,
I map2vec, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SEAICE_MAP_RS2VEC
C | o maps 2 2D-RS-fields to vector and back
C *==========================================================*
C | written by Martin Losch, Oct 2012
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
C === Routine arguments ===
INTEGER n
LOGICAL map2vec
INTEGER myThid
_RS xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RS yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vector (n,nSx,nSy)
#if (defined SEAICE_ALLOW_JFNK) (defined SEAICE_ALLOW_KRYLOV)
C === local variables ===
INTEGER I, J, bi, bj
INTEGER ii, jj, m
CEOP
m = n/2
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef SEAICE_JFNK_MAP_REORDER
ii = 0
IF ( map2vec ) THEN
DO J=1,sNy
jj = 2*sNx*(J-1)
DO I=1,sNx
ii = jj + 2*I
vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
vector(ii, bi,bj) = yfld2d(I,J,bi,bj)
ENDDO
ENDDO
ELSE
DO J=1,sNy
jj = 2*sNx*(J-1)
DO I=1,sNx
ii = jj + 2*I
xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
yfld2d(I,J,bi,bj) = vector(ii, bi,bj)
ENDDO
ENDDO
ENDIF
#else
IF ( map2vec ) THEN
DO J=1,sNy
jj = sNx*(J-1)
DO I=1,sNx
ii = jj + I
vector(ii, bi,bj) = xfld2d(I,J,bi,bj)
vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
ENDDO
ENDDO
ELSE
DO J=1,sNy
jj = sNx*(J-1)
DO I=1,sNx
ii = jj + I
xfld2d(I,J,bi,bj) = vector(ii, bi,bj)
yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* SEAICE_JFNK_MAP_REORDER */
C bi,bj-loops
ENDDO
ENDDO
#endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: SEAICE_FGMRES
C !INTERFACE:
SUBROUTINE SEAICE_FGMRES (
I n,im,rhs,
U sol,i,its,vv,w,wk1,wk2,
I eps,maxits,iout,
U icode,
I myThid )
C-----------------------------------------------------------------------
C mlosch Oct 2012: modified the routine further to be compliant with
C MITgcm standards:
C f90 -> F
C !-comment -> C-comment
C add its to list of arguments
C double precision -> _RL
C implicit none
C
C jfl Dec 1st 2006. We modified the routine so that it is double precison.
C Here are the modifications:
C 1) implicit real (a-h,o-z) becomes implicit real*8 (a-h,o-z)
C 2) real bocomes real*8
C 3) subroutine scopy.f has been changed for dcopy.f
C 4) subroutine saxpy.f has been changed for daxpy.f
C 5) function sdot.f has been changed for ddot.f
C 6) 1e-08 becomes 1d-08
C
C Be careful with the dcopy, daxpy and ddot code...there is a slight
C difference with the single precision versions (scopy, saxpy and sdot).
C In the single precision versions, the array are declared sightly differently.
C It is written for single precision:
C
C modified 12/3/93, array(1) declarations changed to array(*)
C-----------------------------------------------------------------------
implicit none
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
integer myThid
integer i, n, im, its, maxits, iout, icode
_RL rhs(n,nSx,nSy), sol(n,nSx,nSy)
_RL vv(n,im+1,nSx,nSy), w(n,im,nSx,nSy)
_RL wk1(n,nSx,nSy), wk2(n,nSx,nSy), eps
C-----------------------------------------------------------------------
C flexible GMRES routine. This is a version of GMRES which allows a
C a variable preconditioner. Implemented with a reverse communication
C protocole for flexibility -
C DISTRIBUTED VERSION (USES DISTDOT FOR DDOT)
C explicit (exact) residual norms for restarts
C written by Y. Saad, modified by A. Malevsky, version February 1, 1995
C-----------------------------------------------------------------------
C This Is A Reverse Communication Implementation.
C-------------------------------------------------
C USAGE: (see also comments for icode below). FGMRES
C should be put in a loop and the loop should be active for as
C long as icode is not equal to 0. On return fgmres will
C 1) either be requesting the new preconditioned vector applied
C to wk1 in case icode.eq.1 (result should be put in wk2)
C 2) or be requesting the product of A applied to the vector wk1
C in case icode.eq.2 (result should be put in wk2)
C 3) or be terminated in case icode .eq. 0.
C on entry always set icode = 0. So icode should be set back to zero
C upon convergence.
C-----------------------------------------------------------------------
C Here is a typical way of running fgmres:
C
C icode = 0
C 1 continue
C call fgmres (n,im,rhs,sol,i,vv,w,wk1,wk2,eps,maxits,iout,
C & icode,its,mythid)
C
C if (icode .eq. 1) then
C call precon(n, wk1, wk2) <--- user variable preconditioning
C goto 1
C else if (icode .ge. 2) then
C call matvec (n,wk1, wk2) <--- user matrix vector product.
C goto 1
C else
C ----- done ----
C .........
C-----------------------------------------------------------------------
C list of parameters
C-------------------
C
C n == integer. the dimension of the problem
C im == size of Krylov subspace: should not exceed 50 in this
C version (can be reset in code. looking at comment below)
C rhs == vector of length n containing the right hand side
C sol == initial guess on input, approximate solution on output
C vv == work space of size n x (im+1)
C w == work space of length n x im
C wk1,
C wk2, == two work vectors of length n each used for the reverse
C communication protocole. When on return (icode .ne. 1)
C the user should call fgmres again with wk2 = precon * wk1
C and icode untouched. When icode.eq.1 then it means that
C convergence has taken place.
C
C eps == tolerance for stopping criterion. process is stopped
C as soon as ( ||.|| is the euclidean norm):
C || current residual||/||initial residual|| <= eps
C
C maxits== maximum number of iterations allowed
C
C i == internal iteration counter, updated in this routine
C its == current (Krylov) iteration counter, updated in this routine
C
C iout == output unit number number for printing intermediate results
C if (iout .le. 0) no statistics are printed.
C
C icode = integer. indicator for the reverse communication protocole.
C ON ENTRY : icode should be set to icode = 0.
C ON RETURN:
C * icode .eq. 1 value means that fgmres has not finished
C and that it is requesting a preconditioned vector before
C continuing. The user must compute M**(-1) wk1, where M is
C the preconditioing matrix (may vary at each call) and wk1 is
C the vector as provided by fgmres upun return, and put the
C result in wk2. Then fgmres must be called again without
C changing any other argument.
C * icode .eq. 2 value means that fgmres has not finished
C and that it is requesting a matrix vector product before
C continuing. The user must compute A * wk1, where A is the
C coefficient matrix and wk1 is the vector provided by
C upon return. The result of the operation is to be put in
C the vector wk2. Then fgmres must be called again without
C changing any other argument.
C * icode .eq. 0 means that fgmres has finished and sol contains
C the approximate solution.
C comment: typically fgmres must be implemented in a loop
C with fgmres being called as long icode is returned with
C a value .ne. 0.
C-----------------------------------------------------------------------
#if (defined SEAICE_ALLOW_JFNK) (defined SEAICE_ALLOW_KRYLOV)
C local variables -- !jfl modif
integer imax
parameter ( imax = 50 )
_RL hh(4*imax+1,4*imax,MAX_NO_THREADS)
_RL c(4*imax,MAX_NO_THREADS),s(4*imax,MAX_NO_THREADS)
_RL rs(4*imax+1,MAX_NO_THREADS),t,ro
C-------------------------------------------------------------
C arnoldi size should not exceed 50 in this version..
C-------------------------------------------------------------
integer i1Thid(MAX_NO_THREADS)
integer i1, ii, j, jj, k, k1!, n1
integer bi, bj
_RL r0, gam, eps1(MAX_NO_THREADS)
CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
C local common blocks to replace a save statement in the original code
COMMON /SEAICE_FMRES_LOC_I/ i1Thid
COMMON /SEAICE_FMRES_LOC_RL/ hh, c, s, rs, eps1
_RL epsmac
parameter ( epsmac = 1.d-16 )
C
C computed goto
C
if ( im .gt. imax ) stop 'size of krylov space > 50'
goto (100,200,300,11) icode +1
100 continue
its = 0
C-------------------------------------------------------------
C ** outer loop starts here..
C--------------compute initial residual vector --------------
C 10 continue
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j=1,n
wk1(j,bi,bj)=sol(j,bi,bj)
enddo
enddo
enddo
icode = 3
RETURN
11 continue
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j=1,n
vv(j,1,bi,bj) = rhs(j,bi,bj) - wk2(j,bi,bj)
enddo
enddo
enddo
20 continue
call SEAICE_SCALPROD(n, im+1, 1, 1, vv, vv, ro, myThid)
ro = sqrt(ro)
if (ro .eq. 0.0 _d 0) goto 999
t = 1.0 _d 0/ ro
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j=1, n
vv(j,1,bi,bj) = vv(j,1,bi,bj)*t
enddo
enddo
enddo
if (its .eq. 0) eps1(myThid)=eps
C not sure what this is, r0 is never used again
if (its .eq. 0) r0 = ro
if (iout .gt. 0) then
_BEGIN_MASTER( myThid )
write(msgBuf, 199) its, ro
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
C print *,'chau',its, ro !write(iout, 199) its, ro
_END_MASTER( myThid )
endif
C
C initialize 1-st term of rhs of hessenberg system..
C
rs(1,myThid) = ro
i = 0
4 continue
i=i+1
its = its + 1
i1Thid(myThid) = i + 1
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1, n
wk1(k,bi,bj) = vv(k,i,bi,bj)
enddo
enddo
enddo
C
C return
C
icode = 1
RETURN
200 continue
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1, n
w(k,i,bi,bj) = wk2(k,bi,bj)
enddo
enddo
enddo
C
C call matvec operation
C
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
wk1(k,bi,bj)=wk2(k,bi,bj)
enddo
enddo
enddo
C
C return
C
icode = 2
RETURN
300 continue
C
C first call to ope corresponds to intialization goto back to 11.
C
C if (icode .eq. 3) goto 11
i1 = i1Thid(myThid)
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
vv(k,i1,bi,bj)=wk2(k,bi,bj)
enddo
enddo
enddo
C
C modified gram - schmidt...
C
do j=1, i
call SEAICE_SCALPROD(n, im+1, j, i1, vv, vv, t, myThid)
hh(j,i,myThid) = t
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
vv(k,i1,bi,bj) = vv(k,i1,bi,bj) - t*vv(k,j,bi,bj)
enddo
enddo
enddo
enddo
call SEAICE_SCALPROD(n, im+1, i1, i1, vv, vv, t, myThid)
t = sqrt(t)
hh(i1,i,myThid) = t
if (t .ne. 0.0 _d 0) then
t = 1.0 _d 0 / t
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
vv(k,i1,bi,bj) = vv(k,i1,bi,bj)*t
enddo
enddo
enddo
endif
C
C done with modified gram schmidt and arnoldi step.
C now update factorization of hh
C
if (i .ne. 1) then
C
C perfrom previous transformations on i-th column of h
C
do k=2,i
k1 = k-1
t = hh(k1,i,myThid)
hh(k1,i,myThid) = c(k1,myThid)*t+s(k1,myThid)*hh(k,i,myThid)
hh(k,i,myThid) = -s(k1,myThid)*t+c(k1,myThid)*hh(k,i,myThid)
enddo
endif
gam = sqrt(hh(i,i,myThid)**2 + hh(i1,i,myThid)**2)
if (gam .eq. 0.0 _d 0) gam = epsmac
C-----------#determine next plane rotation #-------------------
c(i,myThid) = hh(i, i,myThid)/gam
s(i,myThid) = hh(i1,i,myThid)/gam
C numerically more stable Givens rotation, but the results
C are not better
CML c(i,myThid)=1. _d 0
CML s(i,myThid)=0. _d 0
CML if ( abs(hh(i1(myThid,myThid),i)) .gt. 0.0 _d 0) then
CML if ( abs(hh(i1,i,myThid)) .gt. abs(hh(i,i,myThid)) ) then
CML gam = hh(i,i,myThid)/hh(i1,i,myThid)
CML s(i,myThid) = 1./sqrt(1.+gam*gam)
CML c(i,myThid) = s(i,myThid)*gam
CML else
CML gam = hh(i1,i,myThid)/hh(i,i,myThid)
CML c(i,myThid) = 1./sqrt(1.+gam*gam)
CML s(i,myThid) = c(i,myThid)*gam
CML endif
CML endif
rs(i1,myThid) = -s(i,myThid)*rs(i,myThid)
rs(i ,myThid) = c(i,myThid)*rs(i,myThid)
C
C determine res. norm. and test for convergence
C
hh(i,i,myThid) = c(i,myThid)*hh(i, i,myThid)
& + s(i,myThid)*hh(i1,i,myThid)
ro = abs(rs(i1,myThid))
if (iout .gt. 0) then
_BEGIN_MASTER( myThid )
write(msgBuf, 199) its, ro
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
_END_MASTER( myThid )
endif
CML The original code uses only im as the maximum allowed iterations, but
CML I want no more than maxits iterations, so we change the code here.
CML if (i .lt. im .and. (ro .gt. eps1(myThid))) goto 4
if (i .lt. im .and. its .lt. maxits
& .and. (ro .gt. eps1(myThid))) goto 4
C
C now compute solution. first solve upper triangular system.
C
rs(i,myThid) = rs(i,myThid)/hh(i,i,myThid)
do ii=2,i
k=i-ii+1
k1 = k+1
t=rs(k,myThid)
do j=k1,i
t = t-hh(k,j,myThid)*rs(j,myThid)
enddo
rs(k,myThid) = t/hh(k,k,myThid)
enddo
C
C done with back substitution..
C now form linear combination to get solution
C
do j=1, i
t = rs(j,myThid)
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
sol(k,bi,bj) = sol(k,bi,bj) + t*w(k,j,bi,bj)
enddo
enddo
enddo
enddo
C
C test for return
C
if (ro .le. eps1(myThid) .or. its .ge. maxits) goto 999
C
C else compute residual vector and continue..
C
C goto 10
i1 = i1Thid(myThid)
do j=1,i
jj = i1-j+1
rs(jj-1,myThid) = -s(jj-1,myThid)*rs(jj,myThid)
rs(jj ,myThid) = c(jj-1,myThid)*rs(jj,myThid)
enddo
do j=1,i1
t = rs(j,myThid)
if (j .eq. 1) t = t - 1.0 _d 0
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do k=1,n
vv(k,1,bi,bj) = vv(k,1,bi,bj) + t*vv(k,j,bi,bj)
enddo
enddo
enddo
enddo
C
C restart outer loop.
C
goto 20
999 icode = 0
199 format(' SEAICE_FGMRES: its =', i4, ' res. norm =', d26.16)
C
#endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
RETURN
C-----end-of-fgmres-----------------------------------------------------
C-----------------------------------------------------------------------
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: SEAICE_SCALPROD
C !INTERFACE:
subroutine SEAICE_SCALPROD(n,im,i1,i2,dx,dy,t,myThid)
C forms the dot product of two vectors.
C uses unrolled loops for increments equal to one.
C jack dongarra, linpack, 3/11/78.
C ML: code stolen from BLAS-ddot and adapted for parallel applications
implicit none
#include "SIZE.h"
#include "EEPARAMS.h"
#include "EESUPPORT.h"
#include "SEAICE_SIZE.h"
#include "SEAICE.h"
integer n, im, i1, i2
_RL dx(n,im,nSx,nSy),dy(n,im,nSx,nSy)
_RL t
integer myThid
#if (defined SEAICE_ALLOW_JFNK) (defined SEAICE_ALLOW_KRYLOV)
C local arrays
_RL dtemp(nSx,nSy)
integer i,m,mp1,bi,bj
CEOP
m = mod(n,5)
mp1 = m + 1
t = 0. _d 0
c if( m .eq. 0 ) go to 40
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
dtemp(bi,bj) = 0. _d 0
if ( m .ne. 0 ) then
do i = 1,m
dtemp(bi,bj) = dtemp(bi,bj) + dx(i,i1,bi,bj)*dy(i,i2,bi,bj)
& * scalarProductMetric(i,1,bi,bj)
enddo
endif
if ( n .ge. 5 ) then
c if( n .lt. 5 ) go to 60
c40 mp1 = m + 1
do i = mp1,n,5
dtemp(bi,bj) = dtemp(bi,bj) +
& dx(i, i1,bi,bj)*dy(i, i2,bi,bj)
& * scalarProductMetric(i, 1, bi,bj) +
& dx(i + 1,i1,bi,bj)*dy(i + 1,i2,bi,bj)
& * scalarProductMetric(i + 1,1, bi,bj) +
& dx(i + 2,i1,bi,bj)*dy(i + 2,i2,bi,bj)
& * scalarProductMetric(i + 2,1, bi,bj) +
& dx(i + 3,i1,bi,bj)*dy(i + 3,i2,bi,bj)
& * scalarProductMetric(i + 3,1, bi,bj) +
& dx(i + 4,i1,bi,bj)*dy(i + 4,i2,bi,bj)
& * scalarProductMetric(i + 4,1, bi,bj)
enddo
c60 continue
endif
enddo
enddo
CALL GLOBAL_SUM_TILE_RL( dtemp,t,myThid )
#endif /* SEAICE_ALLOW_JFNK or SEAICE_ALLOW_KRYLOV */
RETURN
END