Sources Esope | Cast3M

prgmr2
C PRGMR2    SOURCE    GOUNAND   22/08/25    21:15:11     11434          
      SUBROUTINE PRGMR2(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,INCX,
     $     ITER,INMV,RESID,KPREC,
     $     RESTRT,ICALRS,IDDOT,IMVEC,IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : PRGMR2
C DESCRIPTION :
C     Préparation de la résolution d'un système linéaire Ax=b
C     par une méthode GMRES(m) (restarted Generalized Minimal
C     Residual) préconditionnée ou non.
C
C     A peut etre symétrique (mais mieux vaut utiliser
C                             le gradient conjugué)
C          ou non symétrique (le résidu b-Ax est toujours DECROISSANT
C                             mais peut éventuellement stagner si
C                             m=RESTRT est trop faible).
C
C     Les 4 préconditionnements disponibles pour cette méthode
C     sont :
C       * Jacobi (diagonal)
C       * D-ILU (Diagonal Incomplete LU factorization)
C       * ILU(0) (Incomplete LU factorization
C                 of level zero ie Choleski incomplet)
C       * MILU(0) (Modified ILU(0)) relaxé
C       5 : préconditionnement ILUT (dual truncation)
C       6 : préconditionnement ILUT2 (une variante du
C           précédent qui remplit mieux la mémoire et
C           fonctionne mieux quelquefois)
C
C     Ces préconditionneurs sont supposés déjà construits (ie
C     factorisés), voir les subroutines MEIDIA, PRDILU, PRILU0,
C     PRMILU et PRILUT.
C
C     Ce sous-programme est en fait une interface aux sous-programmes :
C       gmr, gmrd, gmrdi, gmri0
C     qui sont en Fortran presque pur (pour raison de rapidité)
C     et effectuent la résolution du système linéaire
C     proprement dite.
C
C LANGAGE     : ESOPE
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 (Calcul) : GMR, GMRD, GMRDI, GMRI0
C***********************************************************************
C ENTREES            : MATRIK, MAPREC, KPREC, RESTRT, IMPR
C ENTREES/SORTIES    : INCX, ITER, RESID
C SORTIES            : LRES, IRET
C CODE RETOUR (IRET) : 0 si ok
C                     &lt;0 si problème
C                     >0 si l'algorithme appelé n'a pas convergé
C     MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
C              on pioche dedans les informations nécessaires
C              (différents pointeurs, nb. de ddl...)
C     MAPREC : pointeur sur segment MATRIK de l'include SMMATRIK
C              on va l'utiliser comme préconditionneur.
C     KPREC  : type de préconditionnement à effectuer
C              0 : pas de préconditionnement
C              1 : préconditionnement Jacobi (diagonal)
C              2 : préconditionnement D-ILU
C              3 : préconditionnement ILU(0) (Crout incomplet)
C              4 : préconditionnement MILU(0)
C                     (Crout incomplet modifié)
C              5 : préconditionnement ILUT (dual truncation)
C              6 : préconditionnement ILUT2 (une variante du
C                  précédent qui remplit mieux la mémoire et
C                  fonctionne mieux quelquefois)
C     RESTRT : INTEGER
C              paramètre de redémarrage (restart) (i.e m) pour
C              la méthode GMRES(m).
C              En gros, on 'minimise' le résidu par rapport
C              à ITER&lt;(m=RESTRT) vecteurs orthogonaux 2 à 2.
C              A ITER=multiple(m), on réinitialise le tableau
C              de vecteurs orthogonaux à zéro.
C              Plus m est grand, plus la convergence est bonne
C              et plus l'occupation mémoire et le temps de calcul
C              moyen pour une itération sont élevés.
C              A contrario, si m est trop petit, on peut avoir
C              stagnation du résidu.
C     IMPR   : niveau d'impression (0..4)
C     INCX   : pointeur sur segment IZA de l'include SMMATRIK
C              C'est le 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     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    : v1, 16/06/98, version initiale
C HISTORIQUE : v1, 16/06/98, création
C HISTORIQUE : 09/02/99 : ajout préconditionneur &lt;> MATRIK
C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau gmri0
C (factorisations incomplètes stockées au format MSR (Modified Sparse
C Row) (voir la doc de Sparskit version 2+)
C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
C HISTORIQUE : 09/04/04 : ajout des préconditionneurs ILUTP
C HISTORIQUE : 10/02/06 : Nettoyage
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 et pointeurs associés ..
 
-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
      INTEGER JG
      POINTEUR LRES.MLREEL
-INC SMLENTI
      POINTEUR LNMV.MLENTI
-INC SMMATRIK
      POINTEUR MAPREC.MATRIK
      POINTEUR KMORS.PMORS
      POINTEUR KISA.IZA
      POINTEUR KS2B.IZA
      POINTEUR INCX.IZA
      POINTEUR INVDIA.IZA
      POINTEUR ILUM.PMORS
      POINTEUR ILUI.IZA
*     .. Parameters
*     This is a restart parameter for GMRES
      INTEGER  RESTRT
      INTEGER LDW,NLDW
      SEGMENT SPACE
      REAL*8  WRK(LDW,NLDW)
      ENDSEGMENT
      SEGMENT SPACE2
      POINTEUR WRK2(NIZA).IZA
      ENDSEGMENT
*     .. Work Vectors for GMRES..
      POINTEUR ITMP.IZA
      POINTEUR IV.SPACE2
      POINTEUR IH.SPACE
      POINTEUR IR.IZA,IW.IZA,IAV.IZA
      POINTEUR IY.IZA,IGCOS.IZA,IGSIN.IZA
      POINTEUR IS.IZA
*     .. Scalar Arguments ..
      INTEGER ITER, KPREC, IMPR, IRET
      REAL*8  RESID
      INTEGER NBVA,NJA,NTT,NTTPRE
*     .. Executable Statements ..
*
      IRET = 0
*
* On récupère les paramètres
*
      SEGACT MATRIK
      SEGACT MAPREC
      IF (KSYM.EQ.0) THEN
         IF (IMPR.GT.0) THEN
            WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
     $           'use a Conjugate Gradient instead !'
         ENDIF
C         IRET=-2
C         GOTO 9999
      ENDIF
*     Pour le préconditionneur
      ILUM =MAPREC.KIDMAT(6)
      ILUI =MAPREC.KIDMAT(7)
      IDMAT=MAPREC.KIDMAT(1)
      SEGACT IDMAT
      INVDIA=IDIAG
      SEGDES IDMAT
 
      SEGACT KMORS
      NTT =KMORS.IA(/1)-1
*      NJA =KMORS.JA(/1)
      SEGACT KISA
      SEGACT KS2B
      SEGACT INCX*MOD
C Paramètres des préconditionnements diagonaux et D-ILU
      IF (KPREC.NE.0) THEN
         IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
C     Est-il compatible avec la matrice ?
            SEGACT INVDIA
            NTTPRE=INVDIA.A(/1)
            IF (NTTPRE.NE.NTT) THEN
               WRITE(IOIMP,*) 'La matrice et son préconditionnement'
               WRITE(IOIMP,*) 'ne sont pas compatibles...'
               WRITE(IOIMP,*) 'NTT=',NTT
               WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
               IRET=-2
               GOTO 9999
            ENDIF
C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
C ilutp, ilutpg, ilutpg2
         ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
            SEGACT ILUM
            SEGACT ILUI
            NTTPRE=ILUM.IA(/1)
            IF (NTTPRE.NE.NTT) THEN
               WRITE(IOIMP,*) 'La matrice et son préconditionnement',
     $              'ne sont pas compatibles...'
               WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
               IRET=-2
               GOTO 9999
            ENDIF
         ENDIF
      ENDIF
C
C     Initialisation de l'historique de convergence
C
      JG=ITER+1
      SEGINI LNMV
      SEGINI LRES
C
C
C     Initialisation des vecteurs de travail pour GMRES
C
      NBVA=NTT
      SEGINI IR,IW,IAV
      NBVA=RESTRT
      SEGINI IY,IGCOS,IGSIN
      NBVA=RESTRT+1
      SEGINI IS
      NBVA=NTT
      NIZA=RESTRT+1
      SEGINI IV
      DO I=1,NIZA
         SEGINI ITMP
         IV.WRK2(I)=ITMP
      ENDDO
      LDW=RESTRT+1
      NLDW=RESTRT
      SEGINI IH
C
      CALL GMR2(KMORS,KISA,KS2B,INCX,
     $     KPREC,INVDIA,ILUM,ILUI,
     $     LRES,LNMV,
     $     IR,IW,IAV,
     $     IY,IGCOS,IGSIN,IS,
     $     IV, IH,
     $     ITER,INMV,RESTRT,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
* Gestion du CTRL-C
      if (ierr.NE.0) return
C
C     Désactivation-suppression
C
      SEGSUP IH
      DO I=1,NIZA
         ITMP=IV.WRK2(I)
         SEGSUP ITMP
      ENDDO
      SEGSUP IV
      SEGSUP IS
      SEGSUP IY,IGCOS,IGSIN
      SEGSUP IR,IW,IAV
      JG=ITER+1
      SEGADJ LRES
      SEGDES LRES
      SEGADJ LNMV
      SEGDES LNMV
      IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
         SEGDES INVDIA
      ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
         SEGDES ILUI
         SEGDES ILUM
      ENDIF
      SEGDES INCX
      SEGDES KS2B
      SEGDES KISA
      SEGDES KMORS
      SEGDES MAPREC
      SEGDES MATRIK
C
C     A problem has occured in the GMRES method
C
      IF (IRET.LT.0) GOTO 9999
*
*     Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in prgmr2.eso'
      RETURN
*
*     End of prgmr2
*
      END