/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_fgmres.F
ViewVC logotype

Diff of /MITgcm_contrib/torge/itd/code/seaice_fgmres.F

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

revision 1.3 by torge, Mon Dec 10 22:19:49 2012 UTC revision 1.4 by torge, Wed Mar 27 18:59:52 2013 UTC
# Line 7  C--   File seaice_fgmres.F: seaice fgmre Line 7  C--   File seaice_fgmres.F: seaice fgmre
7  C--   Contents  C--   Contents
8  C--   o SEAICE_FGMRES_DRIVER  C--   o SEAICE_FGMRES_DRIVER
9  C--   o SEAICE_MAP2VEC  C--   o SEAICE_MAP2VEC
10    C--   o SEAICE_MAP_RS2VEC
11  C--   o SEAICE_FGMRES  C--   o SEAICE_FGMRES
12  C--   o SCALPROD  C--   o SEAICE_SCALPROD
13    
14  CBOP  CBOP
15  C     !ROUTINE: SEAICE_FGMRES_DRIVER  C     !ROUTINE: SEAICE_FGMRES_DRIVER
# Line 19  C     !INTERFACE: Line 20  C     !INTERFACE:
20       U     duIce, dvIce,       U     duIce, dvIce,
21       U     iCode,       U     iCode,
22       I     FGMRESeps, iOutFGMRES,       I     FGMRESeps, iOutFGMRES,
23       I     newtonIter, krylovIter, myTime, myIter, myThid )       I     newtonIter,
24         U     krylovIter,
25         I     myTime, myIter, myThid )
26    
27  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
28  C     *==========================================================*  C     *==========================================================*
# Line 46  C     === Routine arguments === Line 49  C     === Routine arguments ===
49  C     myTime :: Simulation time  C     myTime :: Simulation time
50  C     myIter :: Simulation timestep number  C     myIter :: Simulation timestep number
51  C     myThid :: my Thread Id. number  C     myThid :: my Thread Id. number
52  C     newtonIter :: current iterate of Newton iteration  C     newtonIter :: current iterate of Newton iteration (for diagnostics)
53  C     krylovIter :: current iterate of Newton iteration  C     krylovIter :: current iterate of Newton iteration (updated)
54  C     iCode      :: FGMRES parameter to determine next step  C     iCode      :: FGMRES parameter to determine next step
55  C     iOutFGMRES :: control output of fgmres  C     iOutFGMRES :: control output of fgmres
56        _RL     myTime        _RL     myTime
# Line 71  C     u/vIceRes :: residual F(u) Line 74  C     u/vIceRes :: residual F(u)
74        defined (SEAICE_ALLOW_DYNAMICS) )        defined (SEAICE_ALLOW_DYNAMICS) )
75  C     Local variables:  C     Local variables:
76  C     k :: loop indices  C     k :: loop indices
77        INTEGER k        INTEGER k, bi, bj
78  C     FGMRES parameters  C     FGMRES parameters
79  C     n  :: size of the input vector(s)  C     nVec    :: size of the input vector(s)
80  C     im :: size of Krylov space  C     im      :: size of Krylov space
81  C     ifgmres :: interation counter  C     ifgmres :: interation counter
82        INTEGER n        INTEGER nVec
83        PARAMETER ( n  = 2*sNx*sNy*nSx*nSy )        PARAMETER ( nVec  = 2*sNx*sNy )
84        INTEGER im        INTEGER im
85        PARAMETER ( im = 50 )        PARAMETER ( im = 50 )
86        INTEGER ifgmres        INTEGER ifgmres
87  C     work arrays  C     work arrays
88        _RL rhs(n), sol(n)        _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy)
89        _RL vv(n,im+1), w(n,im)        _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy)
90        _RL wk1(n), wk2(n)        _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy)
91  C     need to store some of the fgmres parameters and fields so that  C     need to store some of the fgmres parameters and fields so that
92  C     they are not forgotten between Krylov iterations  C     they are not forgotten between Krylov iterations
93        COMMON /FGMRES_I/ ifgmres        COMMON /FGMRES_I/ ifgmres
94        COMMON /FGMRES_RL/ sol, rhs, vv, w        COMMON /FGMRES_RL/ sol, rhs, vv, w
95  CEOP  CEOP
96    
 C     For now, let only the master thread do all the work  
 C     - copy from 2D arrays to 1D-vector  
 C     - perform fgmres step (including global sums)  
 C     - copy from 1D-vector to 2D arrays  
 C     not sure if this works properly  
   
       _BEGIN_MASTER ( myThid )  
97        IF ( iCode .EQ. 0 ) THEN        IF ( iCode .EQ. 0 ) THEN
98  C     The first guess is zero because it is a correction, but this  C     The first guess is zero because it is a correction, but this
99  C     is implemented by setting du/vIce=0 outside of this routine;  C     is implemented by setting du/vIce=0 outside of this routine;
100  C     this make it possible to restart FGMRES with a nonzero sol  C     this make it possible to restart FGMRES with a nonzero sol
101         CALL SEAICE_MAP2VEC(n,duIce,dvIce,sol,.TRUE.,myThid)         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid)
102  C     wk2 needs to be reset for iCode = 0, because it may contain  C     wk2 needs to be reset for iCode = 0, because it may contain
103  C     remains of the previous Krylov iteration  C     remains of the previous Krylov iteration
104         DO k=1,n         DO bj=myByLo(myThid),myByHi(myThid)
105          wk2(k) = 0. _d 0          DO bi=myBxLo(myThid),myBxHi(myThid)
106             DO k=1,nVec
107              wk2(k,bi,bj) = 0. _d 0
108             ENDDO
109            ENDDO
110         ENDDO         ENDDO
111        ELSEIF ( iCode .EQ. 3 ) THEN        ELSEIF ( iCode .EQ. 3 ) THEN
112         CALL SEAICE_MAP2VEC(n,uIceRes,vIceRes,rhs,.TRUE.,myThid)         CALL SEAICE_MAP2VEC(nVec,uIceRes,vIceRes,rhs,.TRUE.,myThid)
113  C     change sign of rhs because we are solving J*u = -F  C     change sign of rhs because we are solving J*u = -F
114  C     wk2 needs to be initialised for iCode = 3, because it may contain  C     wk2 needs to be initialised for iCode = 3, because it may contain
115  C     garbage  C     garbage
116         DO k=1,n         DO bj=myByLo(myThid),myByHi(myThid)
117          rhs(k) = -rhs(k)          DO bi=myBxLo(myThid),myBxHi(myThid)
118          wk2(k) = 0. _d 0           DO k=1,nVec
119              rhs(k,bi,bj) = -rhs(k,bi,bj)
120              wk2(k,bi,bj) = 0. _d 0
121             ENDDO
122            ENDDO
123         ENDDO         ENDDO
124        ELSE        ELSE
125  C     map preconditioner results or Jacobian times vector,  C     map preconditioner results or Jacobian times vector,
126  C     stored in du/vIce to wk2  C     stored in du/vIce to wk2
127         CALL SEAICE_MAP2VEC(n,duIce,dvIce,wk2,.TRUE.,myThid)         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid)
128        ENDIF        ENDIF
129  C  C
130        CALL SEAICE_FGMRES (n,im,rhs,sol,ifgmres,vv,w,wk1,wk2,        CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,
131       &     FGMRESeps,SEAICEkrylovIterMax,       U     vv,w,wk1,wk2,
132       &     iOutFGMRES,iCode,krylovIter,myThid)       I     FGMRESeps,SEAICEkrylovIterMax,iOutFGMRES,
133         U     iCode,
134         I     myThid)
135  C  C
136        IF ( iCode .EQ. 0 ) THEN        IF ( iCode .EQ. 0 ) THEN
137  C     map sol(ution) vector to du/vIce  C     map sol(ution) vector to du/vIce
138         CALL SEAICE_MAP2VEC(n,duIce,dvIce,sol,.FALSE.,myThid)         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid)
139        ELSE        ELSE
140  C     map work vector to du/vIce to either compute a preconditioner  C     map work vector to du/vIce to either compute a preconditioner
141  C     solution (wk1=rhs) or a Jacobian times wk1  C     solution (wk1=rhs) or a Jacobian times wk1
142         CALL SEAICE_MAP2VEC(n,duIce,dvIce,wk1,.FALSE.,myThid)         CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid)
143        ENDIF        ENDIF
       _END_MASTER ( myThid )  
144    
145  C     Fill overlaps in updated fields  C     Fill overlaps in updated fields
146        CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)        CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid)
# Line 175  C     === Routine arguments === Line 180  C     === Routine arguments ===
180        INTEGER myThid        INTEGER myThid
181        _RL xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
182        _RL yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
183        _RL vector (n)        _RL vector (n,nSx,nSy)
184  C     === local variables ===  C     === local variables ===
185        INTEGER I, J, bi, bj        INTEGER I, J, bi, bj
186        INTEGER ii, jj, ib, jb, m        INTEGER ii, jj, m
187  CEOP  CEOP
188    
189        m = n/2        m = n/2
190        IF ( map2vec ) THEN        DO bj=myByLo(myThid),myByHi(myThid)
191         DO bj=myByLo(myThid),myByHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
192          jb = nSx*sNy*sNx*(bj-1)  #ifdef SEAICE_JFNK_MAP_REORDER
193          DO bi=myBxLo(myThid),myBxHi(myThid)          ii = 0
194           ib = jb + sNy*sNx*(bi-1)          IF ( map2vec ) THEN
195             DO J=1,sNy
196              jj = 2*sNx*(J-1)
197              DO I=1,sNx
198               ii = jj + 2*I
199               vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
200               vector(ii,  bi,bj) = yfld2d(I,J,bi,bj)
201              ENDDO
202             ENDDO
203            ELSE
204           DO J=1,sNy           DO J=1,sNy
205            jj = ib + sNx*(J-1)            jj = 2*sNx*(J-1)
206              DO I=1,sNx
207               ii = jj + 2*I
208               xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
209               yfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
210              ENDDO
211             ENDDO
212            ENDIF
213    #else
214            IF ( map2vec ) THEN
215             DO J=1,sNy
216              jj = sNx*(J-1)
217            DO I=1,sNx            DO I=1,sNx
218             ii = jj + I             ii = jj + I
219             vector(ii)   = xfld2d(I,J,bi,bj)             vector(ii,  bi,bj) = xfld2d(I,J,bi,bj)
220             vector(ii+m) = yfld2d(I,J,bi,bj)             vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
221            ENDDO            ENDDO
222           ENDDO           ENDDO
223          ENDDO          ELSE
224             DO J=1,sNy
225              jj = sNx*(J-1)
226              DO I=1,sNx
227               ii = jj + I
228               xfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
229               yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
230              ENDDO
231             ENDDO
232            ENDIF
233    #endif /* SEAICE_JFNK_MAP_REORDER */
234    C     bi,bj-loops
235         ENDDO         ENDDO
236        ELSE        ENDDO
237         DO bj=myByLo(myThid),myByHi(myThid)  
238          jb = nSx*sNy*sNx*(bj-1)        RETURN
239          DO bi=myBxLo(myThid),myBxHi(myThid)        END
240           ib = jb + sNy*sNx*(bi-1)  
241    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
242    CBOP
243    C     !ROUTINE: SEAICE_MAP_RS2VEC
244    C     !INTERFACE:
245    
246          SUBROUTINE SEAICE_MAP_RS2VEC(
247         I     n,
248         O     xfld2d, yfld2d,
249         U     vector,
250         I     map2vec, myThid )
251    
252    C     !DESCRIPTION: \bv
253    C     *==========================================================*
254    C     | SUBROUTINE SEAICE_MAP_RS2VEC
255    C     | o maps 2 2D-RS-fields to vector and back
256    C     *==========================================================*
257    C     | written by Martin Losch, Oct 2012
258    C     *==========================================================*
259    C     \ev
260    
261    C     !USES:
262          IMPLICIT NONE
263    
264    C     === Global variables ===
265    #include "SIZE.h"
266    #include "EEPARAMS.h"
267    C     === Routine arguments ===
268          INTEGER n
269          LOGICAL map2vec
270          INTEGER myThid
271          _RS xfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
272          _RS yfld2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
273          _RL vector (n,nSx,nSy)
274    C     === local variables ===
275          INTEGER I, J, bi, bj
276          INTEGER ii, jj, m
277    CEOP
278    
279          m = n/2
280          DO bj=myByLo(myThid),myByHi(myThid)
281           DO bi=myBxLo(myThid),myBxHi(myThid)
282    #ifdef SEAICE_JFNK_MAP_REORDER
283            ii = 0
284            IF ( map2vec ) THEN
285           DO J=1,sNy           DO J=1,sNy
286            jj = ib + sNx*(J-1)            jj = 2*sNx*(J-1)
287              DO I=1,sNx
288               ii = jj + 2*I
289               vector(ii-1,bi,bj) = xfld2d(I,J,bi,bj)
290               vector(ii,  bi,bj) = yfld2d(I,J,bi,bj)
291              ENDDO
292             ENDDO
293            ELSE
294             DO J=1,sNy
295              jj = 2*sNx*(J-1)
296              DO I=1,sNx
297               ii = jj + 2*I
298               xfld2d(I,J,bi,bj) = vector(ii-1,bi,bj)
299               yfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
300              ENDDO
301             ENDDO
302            ENDIF
303    #else
304            IF ( map2vec ) THEN
305             DO J=1,sNy
306              jj = sNx*(J-1)
307            DO I=1,sNx            DO I=1,sNx
308             ii = jj + I             ii = jj + I
309             xfld2d(I,J,bi,bj) = vector(ii)             vector(ii,  bi,bj) = xfld2d(I,J,bi,bj)
310             yfld2d(I,J,bi,bj) = vector(ii+m)             vector(ii+m,bi,bj) = yfld2d(I,J,bi,bj)
311            ENDDO            ENDDO
312           ENDDO           ENDDO
313          ENDDO          ELSE
314             DO J=1,sNy
315              jj = sNx*(J-1)
316              DO I=1,sNx
317               ii = jj + I
318               xfld2d(I,J,bi,bj) = vector(ii,  bi,bj)
319               yfld2d(I,J,bi,bj) = vector(ii+m,bi,bj)
320              ENDDO
321             ENDDO
322            ENDIF
323    #endif /* SEAICE_JFNK_MAP_REORDER */
324    C     bi,bj-loops
325         ENDDO         ENDDO
326        ENDIF        ENDDO
327    
328        RETURN        RETURN
329        END        END
# Line 221  C---+----1----+----2----+----3----+----4 Line 332  C---+----1----+----2----+----3----+----4
332  CBOP  CBOP
333  C     !ROUTINE: SEAICE_FGMRES  C     !ROUTINE: SEAICE_FGMRES
334  C     !INTERFACE:  C     !INTERFACE:
335          SUBROUTINE SEAICE_FGMRES (
336        SUBROUTINE SEAICE_FGMRES (n,im,rhs,sol,i,vv,w,wk1, wk2,       I     n,im,rhs,
337       &     eps,maxits,iout,icode,its,myThid)       U     sol,i,its,vv,w,wk1,wk2,
338         I     eps,maxits,iout,
339         U     icode,
340         I     myThid )
341    
342  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
343  C mlosch Oct 2012: modified the routine further to be compliant with  C mlosch Oct 2012: modified the routine further to be compliant with
344  C MITgcm standards:  C MITgcm standards:
345  C f90 -> F  C f90 -> F
346  C !-comment -> C-comment  C !-comment -> C-comment
347    C add its to list of arguments
348  C double precision -> _RL  C double precision -> _RL
349  C implicit none  C implicit none
350  C  C
# Line 251  C modified 12/3/93, array(1) declaration Line 366  C modified 12/3/93, array(1) declaration
366  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
367    
368        implicit none        implicit none
369    C     === Global variables ===
370    #include "SIZE.h"
371    #include "EEPARAMS.h"
372  CML   implicit double precision (a-h,o-z) !jfl modification  CML   implicit double precision (a-h,o-z) !jfl modification
373        integer myThid        integer myThid
374        integer n, im, maxits, iout, icode        integer n, im, its, maxits, iout, icode
375        _RL rhs(*), sol(*), vv(n,im+1), w(n,im)        _RL rhs(n,nSx,nSy), sol(n,nSx,nSy)
376        _RL wk1(n), wk2(n), eps        _RL vv(n,im+1,nSx,nSy), w(n,im,nSx,nSy)
377          _RL wk1(n,nSx,nSy), wk2(n,nSx,nSy), eps
378  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
379  C flexible GMRES routine. This is a version of GMRES which allows a  C flexible GMRES routine. This is a version of GMRES which allows a
380  C a variable preconditioner. Implemented with a reverse communication  C a variable preconditioner. Implemented with a reverse communication
# Line 281  C Here is a typical way of running fgmre Line 400  C Here is a typical way of running fgmre
400  C  C
401  C      icode = 0  C      icode = 0
402  C 1    continue  C 1    continue
403  C      call fgmres (n,im,rhs,sol,i,vv,w,wk1, wk2,eps,maxits,iout,icode)  C      call fgmres (n,im,rhs,sol,i,vv,w,wk1,wk2,eps,maxits,iout,
404    C     &             icode,its,mythid)
405  C  C
406  C      if (icode .eq. 1) then  C      if (icode .eq. 1) then
407  C         call  precon(n, wk1, wk2)    <--- user variable preconditioning  C         call  precon(n, wk1, wk2)    <--- user variable preconditioning
# Line 316  C          || current residual||/||initi Line 436  C          || current residual||/||initi
436  C  C
437  C maxits== maximum number of iterations allowed  C maxits== maximum number of iterations allowed
438  C  C
439    C i     == internal iteration counter, updated in this routine
440    C its   == current (Krylov) iteration counter, updated in this routine
441    C
442  C iout  == output unit number number for printing intermediate results  C iout  == output unit number number for printing intermediate results
443  C          if (iout .le. 0) no statistics are printed.  C          if (iout .le. 0) no statistics are printed.
444  C  C
# Line 350  C     local variables -- !jfl modif Line 473  C     local variables -- !jfl modif
473  C-------------------------------------------------------------  C-------------------------------------------------------------
474  C     arnoldi size should not exceed 50 in this version..  C     arnoldi size should not exceed 50 in this version..
475  C-------------------------------------------------------------  C-------------------------------------------------------------
476        integer i, its, i1, ii, j, jj, k, k1!, n1        integer i, i1, ii, j, jj, k, k1!, n1
477          integer bi, bj
478        _RL r0, gam, epsmac, eps1        _RL r0, gam, epsmac, eps1
479          CHARACTER*(MAX_LEN_MBUF) msgBuf
480    
481  CEOP  CEOP
482        save  CML      save
483    C     local common block to replace the save statement
484          COMMON /SEAICE_FMRES_LOC_I/ i1
485          COMMON /SEAICE_FMRES_LOC_RL/
486         &     hh, c, s, rs, t, ro, r0, gam, epsmac, eps1
487        data epsmac/1.d-16/        data epsmac/1.d-16/
488  C  C
489  C     computed goto  C     computed goto
# Line 362  C Line 491  C
491        if ( im .gt. imax ) stop 'size of krylov space > 50'        if ( im .gt. imax ) stop 'size of krylov space > 50'
492        goto (100,200,300,11) icode +1        goto (100,200,300,11) icode +1
493   100  continue   100  continue
494  CML      n1 = n + 1  CML   n1 = n + 1
495        its = 0        its = 0
496  C-------------------------------------------------------------  C-------------------------------------------------------------
497  C     **  outer loop starts here..  C     **  outer loop starts here..
498  C--------------compute initial residual vector --------------  C--------------compute initial residual vector --------------
499  C 10   continue  C 10   continue
500  CML      call dcopy (n, sol, 1, wk1, 1) !jfl modification  CML   call dcopy (n, sol, 1, wk1, 1) !jfl modification
501        do k=1,n        do bj=myByLo(myThid),myByHi(myThid)
502         wk1(k)=sol(k)         do bi=myBxLo(myThid),myBxHi(myThid)
503            do j=1,n
504             wk1(j,bi,bj)=sol(j,bi,bj)
505            enddo
506           enddo
507        enddo        enddo
508        icode = 3        icode = 3
509        RETURN        RETURN
510   11   continue   11   continue
511        do j=1,n        do bj=myByLo(myThid),myByHi(myThid)
512           vv(j,1) = rhs(j) - wk2(j)         do bi=myBxLo(myThid),myBxHi(myThid)
513            do j=1,n
514             vv(j,1,bi,bj) = rhs(j,bi,bj) - wk2(j,bi,bj)
515            enddo
516           enddo
517        enddo        enddo
518  CML 20   ro = ddot(n, vv, 1, vv,1) !jfl modification   20   continue
519   20   call scalprod(n, vv, vv, ro, myThid)  CML   ro = ddot(n, vv, 1, vv,1) !jfl modification
520          call SEAICE_SCALPROD(n, im+1, 1, 1, vv, vv, ro, myThid)
521        ro = sqrt(ro)        ro = sqrt(ro)
522        if (ro .eq. 0.0d0) goto 999        if (ro .eq. 0.0 _d 0) goto 999
523        t = 1.0d0/ ro        t = 1.0 _d 0/ ro
524        do j=1, n        do bj=myByLo(myThid),myByHi(myThid)
525           vv(j,1) = vv(j,1)*t         do bi=myBxLo(myThid),myBxHi(myThid)
526            do j=1, n
527             vv(j,1,bi,bj) = vv(j,1,bi,bj)*t
528            enddo
529           enddo
530        enddo        enddo
531        if (its .eq. 0) eps1=eps        if (its .eq. 0) eps1=eps
532    C     not sure what this is, r0 is never used again
533        if (its .eq. 0) r0 = ro        if (its .eq. 0) r0 = ro
534        if (iout .gt. 0) write(*, 199) its, ro!&        if (iout .gt. 0) then
535           _BEGIN_MASTER( myThid )
536           write(msgBuf, 199) its, ro
537           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
538         &      SQUEEZE_RIGHT, myThid )
539  C           print *,'chau',its, ro !write(iout, 199) its, ro  C           print *,'chau',its, ro !write(iout, 199) its, ro
540           _END_MASTER( myThid )
541          endif
542  C  C
543  C     initialize 1-st term  of rhs of hessenberg system..  C     initialize 1-st term  of rhs of hessenberg system..
544  C  C
545        rs(1) = ro        rs(1) = ro
546        i = 0        i = 0
547   4    i=i+1   4    continue
548          i=i+1
549        its = its + 1        its = its + 1
550        i1 = i + 1        i1 = i + 1
551        do k=1, n        do bj=myByLo(myThid),myByHi(myThid)
552           wk1(k) = vv(k,i)         do bi=myBxLo(myThid),myBxHi(myThid)
553            do k=1, n
554             wk1(k,bi,bj) = vv(k,i,bi,bj)
555            enddo
556           enddo
557        enddo        enddo
558  C  C
559  C     return  C     return
560  C  C
561        icode = 1        icode = 1
   
562        RETURN        RETURN
563   200  continue   200  continue
564        do k=1, n        do bj=myByLo(myThid),myByHi(myThid)
565           w(k,i) = wk2(k)         do bi=myBxLo(myThid),myBxHi(myThid)
566            do k=1, n
567             w(k,i,bi,bj) = wk2(k,bi,bj)
568            enddo
569           enddo
570        enddo        enddo
571  C  C
572  C     call matvec operation  C     call matvec operation
573  C  C
574        icode = 2  CML   call dcopy(n, wk2, 1, wk1, 1) !jfl modification
575  CML      call dcopy(n, wk2, 1, wk1, 1) !jfl modification        do bj=myByLo(myThid),myByHi(myThid)
576        do k=1,n         do bi=myBxLo(myThid),myBxHi(myThid)
577         wk1(k)=wk2(k)          do k=1,n
578             wk1(k,bi,bj)=wk2(k,bi,bj)
579            enddo
580           enddo
581        enddo        enddo
582  C  C
583  C     return  C     return
584  C  C
585          icode = 2
586        RETURN        RETURN
587   300  continue   300  continue
588  C  C
589  C     first call to ope corresponds to intialization goto back to 11.  C     first call to ope corresponds to intialization goto back to 11.
590  C  C
591  C      if (icode .eq. 3) goto 11  C      if (icode .eq. 3) goto 11
592  CML      call  dcopy (n, wk2, 1, vv(1,i1), 1) !jfl modification  CML   call  dcopy (n, wk2, 1, vv(1,i1), 1) !jfl modification
593        do k=1,n        do bj=myByLo(myThid),myByHi(myThid)
594         vv(k,i1)=wk2(k)         do bi=myBxLo(myThid),myBxHi(myThid)
595            do k=1,n
596             vv(k,i1,bi,bj)=wk2(k,bi,bj)
597            enddo
598           enddo
599        enddo        enddo
600  C  C
601  C     modified gram - schmidt...  C     modified gram - schmidt...
602  C  C
603        do j=1, i        do j=1, i
604  CML         t = ddot(n, vv(1,j), 1, vv(1,i1), 1) !jfl modification  CML    t = ddot(n, vv(1,j), 1, vv(1,i1), 1) !jfl modification
605           call scalprod(n, vv(1,j), vv(1,i1), t, myThid)         call SEAICE_SCALPROD(n, im+1, j, i1, vv, vv, t, myThid)
606           hh(j,i) = t         hh(j,i) = t
607  CML         call daxpy(n, -t, vv(1,j), 1, vv(1,i1), 1) !jfl modification  CML    call daxpy(n, -t, vv(1,j), 1, vv(1,i1), 1) !jfl modification
608  CML      enddo  CML   enddo
609  CML      do j=1, i  CML   do j=1, i
610  CML         t = hh(j,i)  CML    t = hh(j,i)
611           do bj=myByLo(myThid),myByHi(myThid)
612            do bi=myBxLo(myThid),myBxHi(myThid)
613           do k=1,n           do k=1,n
614            vv(k,i1) = vv(k,i1) - t*vv(k,j)            vv(k,i1,bi,bj) = vv(k,i1,bi,bj) - t*vv(k,j,bi,bj)
615           enddo           enddo
616            enddo
617           enddo
618        enddo        enddo
619  CML      t = sqrt(ddot(n, vv(1,i1), 1, vv(1,i1), 1)) !jfl modification  CML   t = sqrt(ddot(n, vv(1,i1), 1, vv(1,i1), 1)) !jfl modification
620        call scalprod(n, vv(1,i1), vv(1,i1), t, myThid)        call SEAICE_SCALPROD(n, im+1, i1, i1, vv, vv, t, myThid)
621        t = sqrt(t)        t = sqrt(t)
622        hh(i1,i) = t        hh(i1,i) = t
623        if (t .eq. 0.0d0) goto 58        if (t .ne. 0.0 _d 0) then
624        t = 1.0d0 / t         t = 1.0 _d 0 / t
625        do k=1,n         do bj=myByLo(myThid),myByHi(myThid)
626           vv(k,i1) = vv(k,i1)*t          do bi=myBxLo(myThid),myBxHi(myThid)
627        enddo           do k=1,n
628              vv(k,i1,bi,bj) = vv(k,i1,bi,bj)*t
629             enddo
630            enddo
631           enddo
632          endif
633  C  C
634  C     done with modified gram schimd and arnoldi step.  C     done with modified gram schimd and arnoldi step.
635  C     now  update factorization of hh  C     now  update factorization of hh
636  C  C
637   58   if (i .eq. 1) goto 121        if (i .ne. 1) then
638  C  C
639  C     perfrom previous transformations  on i-th column of h  C     perfrom previous transformations  on i-th column of h
640  C  C
641        do k=2,i         do k=2,i
642           k1 = k-1          k1 = k-1
643           t = hh(k1,i)          t = hh(k1,i)
644           hh(k1,i) = c(k1)*t + s(k1)*hh(k,i)          hh(k1,i) = c(k1)*t + s(k1)*hh(k,i)
645           hh(k,i) = -s(k1)*t + c(k1)*hh(k,i)          hh(k,i) = -s(k1)*t + c(k1)*hh(k,i)
646        enddo         enddo
647   121  gam = sqrt(hh(i,i)**2 + hh(i1,i)**2)        endif
648        if (gam .eq. 0.0d0) gam = epsmac        gam = sqrt(hh(i,i)**2 + hh(i1,i)**2)
649          if (gam .eq. 0.0 _d 0) gam = epsmac
650  C-----------#determine next plane rotation  #-------------------  C-----------#determine next plane rotation  #-------------------
651        c(i) = hh(i,i)/gam        c(i)   = hh(i,i)/gam
652        s(i) = hh(i1,i)/gam        s(i)   = hh(i1,i)/gam
653    C     numerically more stable Givens rotation, but the results
654    C     are not better
655    CML      c(i)=1. _d 0
656    CML      s(i)=0. _d 0
657    CML      if ( abs(hh(i1,i)) .gt. 0.0 _d 0) then
658    CML       if ( abs(hh(i1,i)) .gt. abs(hh(i,i)) ) then
659    CML        gam = hh(i,i)/hh(i1,i)
660    CML        s(i) = 1./sqrt(1.+gam*gam)
661    CML        c(i) = s(i)*gam
662    CML       else
663    CML        gam = hh(i1,i)/hh(i,i)
664    CML        c(i) = 1./sqrt(1.+gam*gam)
665    CML        s(i) = c(i)*gam
666    CML       endif
667    CML      endif
668        rs(i1) = -s(i)*rs(i)        rs(i1) = -s(i)*rs(i)
669        rs(i) =  c(i)*rs(i)        rs(i)  =  c(i)*rs(i)
670  C  C
671  C     determine res. norm. and test for convergence-  C     determine res. norm. and test for convergence
672  C  C
673        hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i)        hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i)
674        ro = abs(rs(i1))        ro = abs(rs(i1))
675        if (iout .gt. 0) write(*, 199) its, ro        if (iout .gt. 0) then
676           _BEGIN_MASTER( myThid )
677           write(msgBuf, 199) its, ro
678           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
679         &      SQUEEZE_RIGHT, myThid )
680           _END_MASTER( myThid )
681          endif
682        if (i .lt. im .and. (ro .gt. eps1))  goto 4        if (i .lt. im .and. (ro .gt. eps1))  goto 4
683  C  C
684  C     now compute solution. first solve upper triangular system.  C     now compute solution. first solve upper triangular system.
685  C  C
686        rs(i) = rs(i)/hh(i,i)        rs(i) = rs(i)/hh(i,i)
687        do ii=2,i        do ii=2,i
688           k=i-ii+1         k=i-ii+1
689           k1 = k+1         k1 = k+1
690           t=rs(k)         t=rs(k)
691           do j=k1,i         do j=k1,i
692              t = t-hh(k,j)*rs(j)          t = t-hh(k,j)*rs(j)
693           enddo         enddo
694           rs(k) = t/hh(k,k)         rs(k) = t/hh(k,k)
695        enddo        enddo
696  C  C
697  C     done with back substitution..  C     done with back substitution..
# Line 503  C     now form linear combination to get Line 699  C     now form linear combination to get
699  C  C
700        do j=1, i        do j=1, i
701         t = rs(j)         t = rs(j)
702  C         call daxpy(n, t, w(1,j), 1, sol,1) !jfl modification  CML    call daxpy(n, t, w(1,j), 1, sol,1) !jfl modification
703         do k=1,n         do bj=myByLo(myThid),myByHi(myThid)
704          sol(k) = sol(k) + t*w(k,j)          do bi=myBxLo(myThid),myBxHi(myThid)
705             do k=1,n
706              sol(k,bi,bj) = sol(k,bi,bj) + t*w(k,j,bi,bj)
707             enddo
708            enddo
709         enddo         enddo
710        enddo        enddo
711  C  C
# Line 518  C Line 718  C
718  C       goto 10  C       goto 10
719    
720        do j=1,i        do j=1,i
721          jj = i1-j+1         jj = i1-j+1
722          rs(jj-1) = -s(jj-1)*rs(jj)         rs(jj-1) = -s(jj-1)*rs(jj)
723          rs(jj) = c(jj-1)*rs(jj)         rs(jj) = c(jj-1)*rs(jj)
724        enddo        enddo
725        do j=1,i1        do j=1,i1
726          t = rs(j)         t = rs(j)
727          if (j .eq. 1)  t = t-1.0d0         if (j .eq. 1)  t = t-1.0 _d 0
728  CML        call daxpy (n, t, vv(1,j), 1,  vv, 1)  CML    call daxpy (n, t, vv(1,j), 1,  vv, 1)
729          do k=1,n         do bj=myByLo(myThid),myByHi(myThid)
730           vv(k,1) = vv(k,1) + t*vv(k,j)          do bi=myBxLo(myThid),myBxHi(myThid)
731             do k=1,n
732              vv(k,1,bi,bj) = vv(k,1,bi,bj) + t*vv(k,j,bi,bj)
733             enddo
734          enddo          enddo
735           enddo
736        enddo        enddo
737  C  C
738  C     restart outer loop.  C     restart outer loop.
# Line 536  C Line 740  C
740        goto 20        goto 20
741   999  icode = 0   999  icode = 0
742    
743   199  format('   -- fgmres its =', i4, ' res. norm =', d26.16)   199  format(' SEAICE_FGMRES: its =', i4, ' res. norm =', d26.16)
744  C  C
745        RETURN        RETURN
746  C-----end-of-fgmres-----------------------------------------------------  C-----end-of-fgmres-----------------------------------------------------
# Line 545  C--------------------------------------- Line 749  C---------------------------------------
749    
750  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
751  CBOP  CBOP
752  C     !ROUTINE: SCALPROD  C     !ROUTINE: SEAICE_SCALPROD
753  C     !INTERFACE:  C     !INTERFACE:
754    
755        subroutine scalprod(n,dx,dy,t,myThid)        subroutine SEAICE_SCALPROD(n,im,i1,i2,dx,dy,t,myThid)
756    
757  C     forms the dot product of two vectors.  C     forms the dot product of two vectors.
758  C     uses unrolled loops for increments equal to one.  C     uses unrolled loops for increments equal to one.
759  C     jack dongarra, linpack, 3/11/78.  C     jack dongarra, linpack, 3/11/78.
760  C     ML: code stolen from BLAS and adapted for parallel applications  C     ML: code stolen from BLAS-ddot and adapted for parallel applications
761    
762        implicit none        implicit none
763  #include "SIZE.h"  #include "SIZE.h"
764  #include "EEPARAMS.h"  #include "EEPARAMS.h"
765  #include "EESUPPORT.h"  #include "EESUPPORT.h"
766        integer n  #include "SEAICE_SIZE.h"
767        _RL dx(n),dy(n)  #include "SEAICE.h"
768        real*8 t        integer n, im, i1, i2
769          _RL dx(n,im,nSx,nSy),dy(n,im,nSx,nSy)
770          _RL t
771        integer myThid        integer myThid
772    C     local arrays
773        real*8 dtemp        _RL dtemp(nSx,nSy)
774        integer i,m,mp1        integer i,m,mp1,bi,bj
 #ifdef ALLOW_USE_MPI  
       INTEGER mpiRC  
 #endif   /* ALLOW_USE_MPI */  
775  CEOP  CEOP
776  C  
777        m = mod(n,5)        m = mod(n,5)
778        dtemp = 0. _d 0        mp1 = m + 1
779        t     = 0. _d 0        t     = 0. _d 0
780        if( m .eq. 0 ) go to 40  c     if( m .eq. 0 ) go to 40
781        do i = 1,m        do bj=myByLo(myThid),myByHi(myThid)
782         dtemp = dtemp + dx(i)*dy(i)         do bi=myBxLo(myThid),myBxHi(myThid)
783        enddo          dtemp(bi,bj) = 0. _d 0
784        if( n .lt. 5 ) go to 60          if ( m .ne. 0 ) then
785     40 mp1 = m + 1           do i = 1,m
786        do i = mp1,n,5            dtemp(bi,bj) = dtemp(bi,bj) + dx(i,i1,bi,bj)*dy(i,i2,bi,bj)
787         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +       &         * scalarProductMetric(i,1,bi,bj)
788       &      dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) +           enddo
789       &      dx(i + 4)*dy(i + 4)          endif
790        enddo          if ( n .ge. 5 ) then
791     60 continue  c     if( n .lt. 5 ) go to 60
792  C     sum over all processors  c40   mp1 = m + 1
793  #ifdef ALLOW_USE_MPI           do i = mp1,n,5
794        t = dtemp            dtemp(bi,bj) = dtemp(bi,bj)               +
795        IF ( usingMPI ) THEN       &        dx(i,    i1,bi,bj)*dy(i,    i2,bi,bj)
796         CALL MPI_Allreduce(t,dtemp,1,MPI_DOUBLE_PRECISION,MPI_SUM,       &        * scalarProductMetric(i,    1, bi,bj) +
797       &                   MPI_COMM_MODEL,mpiRC)       &        dx(i + 1,i1,bi,bj)*dy(i + 1,i2,bi,bj)
798        ENDIF       &        * scalarProductMetric(i + 1,1, bi,bj) +
799  #endif /* ALLOW_USE_MPI */       &        dx(i + 2,i1,bi,bj)*dy(i + 2,i2,bi,bj)
800        t = dtemp       &        * scalarProductMetric(i + 2,1, bi,bj) +
801         &        dx(i + 3,i1,bi,bj)*dy(i + 3,i2,bi,bj)
802  CML      return       &        * scalarProductMetric(i + 3,1, bi,bj) +
803  CML      end       &        dx(i + 4,i1,bi,bj)*dy(i + 4,i2,bi,bj)
804  CML       &        * scalarProductMetric(i + 4,1, bi,bj)
805  CML      subroutine daxpy(n,da,dx,incx,dy,incy)           enddo
806  CMLC  c60   continue
807  CMLC     constant times a vector plus a vector.          endif
808  CMLC     uses unrolled loops for increments equal to one.         enddo
809  CMLC     jack dongarra, linpack, 3/11/78.        enddo
810  CMLC        CALL GLOBAL_SUM_TILE_RL( dtemp,t,myThid )
 CML      _RL dx(n),dy(n),da  
 CML      integer i,incx,incy,ix,iy,m,mp1,n  
 CMLC  
 CML      if(n.le.0)return  
 CML      if (da .eq. 0.0d0) return  
 CML      if(incx.eq.1.and.incy.eq.1)go to 20  
 CMLC  
 CMLC        code for unequal increments or equal increments  
 CMLC          not equal to 1  
 CMLC  
 CML      ix = 1  
 CML      iy = 1  
 CML      if(incx.lt.0)ix = (-n+1)*incx + 1  
 CML      if(incy.lt.0)iy = (-n+1)*incy + 1  
 CML      do 10 i = 1,n  
 CML        dy(iy) = dy(iy) + da*dx(ix)  
 CML        ix = ix + incx  
 CML        iy = iy + incy  
 CML   10 continue  
 CML      return  
 CMLC  
 CMLC        code for both increments equal to 1  
 CMLC  
 CMLC  
 CMLC        clean-up loop  
 CMLC  
 CML   20 m = mod(n,4)  
 CML      if( m .eq. 0 ) go to 40  
 CML      do 30 i = 1,m  
 CML        dy(i) = dy(i) + da*dx(i)  
 CML   30 continue  
 CML      if( n .lt. 4 ) return  
 CML   40 mp1 = m + 1  
 CML      do 50 i = mp1,n,4  
 CML        dy(i) = dy(i) + da*dx(i)  
 CML        dy(i + 1) = dy(i + 1) + da*dx(i + 1)  
 CML        dy(i + 2) = dy(i + 2) + da*dx(i + 2)  
 CML        dy(i + 3) = dy(i + 3) + da*dx(i + 3)  
 CML   50 continue  
 CML      return  
 CML      end  
 CML  
 CML      subroutine  dcopy(n,dx,incx,dy,incy)  
 CMLC  
 CMLC     copies a vector, x, to a vector, y.  
 CMLC     uses unrolled loops for increments equal to one.  
 CMLC     jack dongarra, linpack, 3/11/78.  
 CMLC  
 CML      _RL dx(n),dy(n)  
 CML      integer i,incx,incy,ix,iy,m,mp1,n  
 CMLC  
 CML      if(n.le.0)return  
 CML      if(incx.eq.1.and.incy.eq.1)go to 20  
 CMLC  
 CMLC        code for unequal increments or equal increments  
 CMLC          not equal to 1  
 CMLC  
 CML      ix = 1  
 CML      iy = 1  
 CML      if(incx.lt.0)ix = (-n+1)*incx + 1  
 CML      if(incy.lt.0)iy = (-n+1)*incy + 1  
 CML      do 10 i = 1,n  
 CML        dy(iy) = dx(ix)  
 CML        ix = ix + incx  
 CML        iy = iy + incy  
 CML   10 continue  
 CML      return  
 CMLC  
 CMLC        code for both increments equal to 1  
 CMLC  
 CMLC  
 CMLC        clean-up loop  
 CMLC  
 CML   20 m = mod(n,7)  
 CML      if( m .eq. 0 ) go to 40  
 CML      do 30 i = 1,m  
 CML        dy(i) = dx(i)  
 CML   30 continue  
 CML      if( n .lt. 7 ) return  
 CML   40 mp1 = m + 1  
 CML      do 50 i = mp1,n,7  
 CML        dy(i) = dx(i)  
 CML        dy(i + 1) = dx(i + 1)  
 CML        dy(i + 2) = dx(i + 2)  
 CML        dy(i + 3) = dx(i + 3)  
 CML        dy(i + 4) = dx(i + 4)  
 CML        dy(i + 5) = dx(i + 5)  
 CML        dy(i + 6) = dx(i + 6)  
 CML   50 continue  
811    
812  #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */  #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */
813    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22