Sources Esope | Cast3M

gmr2
C GMR2      SOURCE    GOUNAND   22/08/25    21:15:04     11434          
      SUBROUTINE GMR2(KMORS,KISA,B,X,
     $     KPREC,INVDIA,ILUM,ILUI,
     $     LRES,LNMV,
     $     R,W,AV,
     $     Y,GCOS,GSIN,S,
     $     V, H,
     $     ITER,INMV,RESTRT,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : GMR2
C DESCRIPTION :
C DESCRIPTION :
C     Résolution d'un système linéaire Ax=b
C     par une méthode GMRES(m) (restarted Generalized Minimal
C     Residual) préconditionné par une factorisation  ILU(0)
C     (Crout incomplet) ou MILU(0) (Crout inc. modifié relaxé) de A.
C
C LANGAGE     : FORTRAN 77 + chouhia ESOPE (pour les E/S)
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
C               mél : gounand@semt2.smts.cea.fr
C REFERENCE (bibtex-like) :
C  @BOOK{templates,
C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
C                H. Van der Vorst},
C        TITLE={Templates for the Solution of Linear Systems :
C               Building Blocks for Iterative Methods},
C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
C  -> URL : http://www.netlib.org/templates/Templates.html
C***********************************************************************
C APPELES (BLAS 1) : DCOPY, DAXPY      (subroutines)
C                    DNRM2, DDOT, DSCAL (functions)
C APPELES (BLAS 2) : DTRSV             (subroutines)
C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
C                    PSILU0 : Preconditionner solve
C                             (Montée-Descente adaptée à ILU(0))
C***********************************************************************
C ENTREES            : N, NNZ,
C                      ROWPTR, COLIND, VAL,
C                      B, JU, JLU, ALU, RESTRT, IMPR
C ENTREES/SORTIES    : X,
C                      R,W,AV,
C                      Y,GCOS,GSIN,S,
C                      V,H,
C                      ITER, RESID
C SORTIES            : LRES, IRET
C CODE RETOUR (IRET) : 0 si ok
C                     &lt;0 si problème (cf. BRTOL)
C                     >0 si l'algorithme n'a pas convergé
C     N      : nombre de degrés de liberté
C     NNZ    : nombre de valeurs non nulles de la matrice Morse
C     ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
C                           et valeurs de la matrice Morse
C     B      : vecteur second membre
C     ILUVAL : valeurs de la matrice de la factorisation
C              ILU(0) ou MILU(0) de A.
C              Cette matrice est stockée en Morse, elle a le
C              meme profil que A.
C     RESTRT : INTEGER
C              paramètre de redémarrage (restart) (i.e. m) pour
C              la méthode GMRES(m). (cf. prgmr.eso)
C     IMPR   : niveau d'impression (0..4)
C     X      : vecteur des inconnues.
C              En entrée, contient l'estimation x(0).
C              En sortie, la solution trouvée (si la méthode
C              a convergé).
C     R,W,AV        : vecteurs de travail initialisés et détruits
C     Y,GCOS,GSIN,S   dans prgmr.
C     V,H           : tableaux de travail...
C     ITER   : type INTEGER.
C              En entrée, contient le nombre maximum
C                d'itérations à effectuer.
C              En sortie, contient le nombre d'itérations
C                réellement effectuées.
C     RESID  : type REAL*8.
C              En entrée, la valeur maximale du résidu à
C                convergence de l'algorithme mesurée comme suit :
C                    norme[L2](b - A*x) / norme[L2]( b )
C              En sortie, la valeur effective de cette mesure.
C     LRES   : pointeur sur segment MLREEL (include SMLREEL)
C              contient l'historique de la valeur de RESID en
C              fonction du nombre d'itérations effectuées.
C
C***********************************************************************
C VERSION    :     20/12/99
C HISTORIQUE : v1, 16/06/98, création
C HISTORIQUE :     20/12/99
C Interfaçage avec le nouveau psilu0 (factorisations incomplètes
C stockées au format MSR (Modified Sparse Row).
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
*
*     .. Includes .. a commenter si mise au point ok
-INC CCREEL
 
-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
      POINTEUR LRES.MLREEL
-INC SMLENTI
      POINTEUR LNMV.MLENTI
-INC SMMATRIK
      POINTEUR KMORS.PMORS
      POINTEUR KISA.IZA
      POINTEUR INVDIA.IZA
      POINTEUR ILUM.PMORS
      POINTEUR ILUI.IZA
      POINTEUR B.IZA
      POINTEUR X.IZA
*     .. Parameters
*     This is a restart parameter for GMRES
      SEGMENT SPACE
      REAL*8  WRK(LDW,NLDW)
      ENDSEGMENT
      SEGMENT SPACE2
      POINTEUR WRK2(NIZA).IZA
      ENDSEGMENT
*     .. Work Vectors for GMRES..
      POINTEUR VTMP.IZA
      POINTEUR V.SPACE2
      POINTEUR H.SPACE
      POINTEUR R.IZA,W.IZA,AV.IZA
      POINTEUR Y.IZA,GCOS.IZA,GSIN.IZA
      POINTEUR S.IZA
 
      INTEGER IMPR
*     .. Scalar Arguments ..
      INTEGER  N,NNZ
      INTEGER  ITER,RESTRT,IRET
      REAL*8   RESID
*     .. Tableaux de travail à supprimer
C      REAL*8  R(N),W(N),AV(N)
C      REAL*8  Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT)
C      REAL*8  S(RESTRT+1)
C      REAL*8  V(N,RESTRT+1)
C      REAL*8  H(RESTRT+1,RESTRT)
*
*
*     .. Variables locales
*     .. Parameters
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D0, ONE=1.0D0)
*     ..
      INTEGER I,MAXIT,J,K
      REAL*8  BNRM2,TOL,GDOT,GNRM2
      REAL*8  H1,H2,RAY
*     ..
*     .. External subroutines and functions..
      EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL,DTRSV
*     ..
*     .. Executable Statements ..
*
*
      IRET  = 0
      INMV  = 0
      MAXIT = ITER
      TOL   = RESID
      IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmr2.eso'
*
*     Set initial residual.
*
C                    r=b
      CALL GCOPY(B,R)
*      IF (GNRM2(X).NE.ZERO) THEN
C                    r=b-Ax
*         CALL GMOMV ('N',ONE,KMORS,KISA,X,ONE,R)
*      ENDIF
      CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
      INMV=INMV+1
*
      ITER = 0
      BNRM2 = GNRM2(B)
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) '||B||=',BNRM2
      ENDIF
      RNRM2 = GNRM2(R)
      IF (ICALRS.EQ.1) BNRM2=RNRM2
*
      IF (BNRM2.LT.XPETIT) BNRM2 = ONE
      RESID = RNRM2 / BNRM2
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'Résidu initial=',RESID
      ENDIF
      LRES.PROG(ITER+1)=RESID
      LNMV.LECT(ITER+1)=INMV
      IF (RESID.LE.TOL) GOTO 30
*
   10 CONTINUE
* Gestion du CTRL-C
      if (ierr.NE.0) return
*
*        Construct the first column of V, and initialize S to the
*        elementary vector E1 scaled by RNORM.
*
C        Preconditionner solve (ILU(0))  : Mv(1)=r
      VTMP=V.WRK2(1)
      CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
     $     ,VTMP,R)
 
C                    s=||v(1)|| e(1)
C                    v(1)=v(1) / ||v(1)||
      S.A(1) = GNRM2(VTMP)
      CALL GSCAL(ONE/S.A(1),VTMP)
*
      DO 50 I = 1, RESTRT
         ITER = ITER + 1
C                    av=A v(i)
         VTMP=V.WRK2(I)
         CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,VTMP,ZERO,AV)
         INMV=INMV+1
C        Preconditionner solve (ILU(0)) : Mw=av
         CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
     $        ,W,AV)
*
*           Construct I-th column of H so that it is orthonormal to
*           the previous I-1 columns.
*
         DO 60 K = 1, I
            VTMP=V.WRK2(K)
            H.WRK(K,I) = GDOT(W,VTMP,IDDOT)
            CALL GAXPY(-ONE*H.WRK(K,I),VTMP,W)
 60      CONTINUE
         H.WRK(I+1,I)=GNRM2(W)
         VTMP=V.WRK2(I+1)
         CALL GCOPY(W,VTMP)
 
         Xval = H.WRK(I+1,I)
         IF(ABS(Xval) .LT. XPETIT )THEN
           CALL GSCAL(XGRAND,VTMP)
         ELSE
           CALL GSCAL(ONE/Xval,VTMP)
         ENDIF
*
C         IF (I.GT.1)
*
*              Apply Givens rotations to the I-th column of H. This
*              effectively reduces the Hessenberg matrix to upper
*              triangular form during the RESTRT iterations.
*
         DO 70 K = 1, I-1
            H1= (GCOS.A(K)*H.WRK(K,I)) + (GSIN.A(K)*H.WRK(K+1,I))
            H2=-(GSIN.A(K)*H.WRK(K,I)) + (GCOS.A(K)*H.WRK(K+1,I))
            H.WRK(K,I)  =H1
            H.WRK(K+1,I)=H2
 70      CONTINUE
 
*
*           Approximate residual norm. Check tolerance. If okay, compute
*           final approximation vector X and quit.
*
         RAY=SQRT((H.WRK(I,I)*H.WRK(I,I)) + (H.WRK(I+1,I)*H.WRK(I+1,I)))
         GCOS.A(I) = H.WRK(I,I) / RAY
         GSIN.A(I) = H.WRK(I+1,I) / RAY
         H.WRK(I,I)=RAY
         H.WRK(I+1,I)=ZERO
         S.A(I+1)=-GSIN.A(I)*S.A(I)
         S.A(I)  = GCOS.A(I)*S.A(I)
         RESID=ABS(S.A(I+1)) / BNRM2
         LRES.PROG(ITER+1)=RESID
         LNMV.LECT(ITER+1)=INMV
*         IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN
         IF (RESID.LE.TOL.OR.INMV.GE.MAXIT) THEN
*                  Solve the system Hy=S
            CALL DCOPY(I,S.A,1,Y.A,1)
            CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H.WRK,RESTRT+1,
     $           Y.A,1)
*                  Updating solution
C            CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
C     $           WORK(1,S),WORK(1,V))
            DO 503 J=I,1,-1
               VTMP=V.WRK2(J)
               CALL GAXPY(Y.A(J),VTMP,X)
 503        CONTINUE
            IF (RESID.LE.TOL) THEN
               GOTO 30
            ENDIF
*            IF (ITER.EQ.MAXIT) THEN
            IF (INMV.GE.MAXIT) THEN
               GOTO 40
            ENDIF
         ENDIF
 50   CONTINUE
*
*        Compute current solution vector X.
*
*                  Solve the system Hy=S
      CALL DCOPY(RESTRT,S.A,1,Y.A,1)
      CALL DTRSV('UPPER','NOTRANS','NONUNIT',
     $     RESTRT,H.WRK,RESTRT+1,Y.A,1)
*                  Updating solution
      DO 303 J=RESTRT,1,-1
         VTMP=V.WRK2(J)
         CALL GAXPY(Y.A(J),VTMP,X)
 303  CONTINUE
*
*        Compute residual vector R, find norm,
*        then check for tolerance.
*
      CALL GCOPY(B,R)
      CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
      INMV=INMV+1
      S.A(RESTRT+1) = GNRM2(R)
      RESID = S.A(RESTRT+1) / BNRM2
      LRES.PROG(ITER+1)=RESID
      LNMV.LECT(ITER+1)=INMV
      IF (RESID.LE.TOL) GOTO 30
*
*     Iteration fails.
*
*      IF (ITER.EQ.MAXIT) THEN
      IF (INMV.GE.MAXIT) THEN
         GOTO 40
      ENDIF
*
*     Do some more
*
      GOTO 10
*
*     Iteration successful
*
   30 CONTINUE
      IRET = 0
      RETURN
*
*     Iteration fails
*
   40 CONTINUE
      IRET = 1
      RETURN
*
*     End of GMR2.
*
      END