Sources Esope | Cast3M

prgmr2
C PRGMR2    SOURCE    CB215821  25/06/19    21:15:02     12288          
      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
C      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
C      SEGDES LRES
      SEGADJ LNMV
C      SEGDES LNMV
C      IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
C         SEGDES INVDIA
C      ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
C         SEGDES ILUI
C         SEGDES ILUM
C      ENDIF
C      SEGDES INCX
C      SEGDES KS2B
C      SEGDES KISA
C      SEGDES KMORS
C      SEGDES MAPREC
C      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