Sources Esope | Cast3M

prcg2
C PRCG2     SOURCE    GOUNAND   22/08/25    21:15:09     11434          
      SUBROUTINE PRCG2(KMORS,KISA,KS2B,MATRIK,MAPREC,
     $     LRES,LNMV,INCX,ITER,INMV,RESID,KPREC,BRTOL,ICALRS,IDDOT,
     $     IMVEC,
     $     IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : PRCG2
C DESCRIPTION :
C     Préparation de la résolution d'un système linéaire Ax=b
C     par la méthode du gradient conjugué (préconditionnée ou non)
C
C     A doit etre symétrique.
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       cg, cgd, cgdilu, cgilu0
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) : CG, CGD, CGDILU, CGILU0
C***********************************************************************
C ENTREES            : MATRIK, MAPREC, KPREC, 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) (Choleski incomplet)
C              4 : préconditionnement MILU(0)
C                     (Choleski 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     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, 02/04/98, version initiale
C HISTORIQUE : v1, 02/04/98, création
C HISTORIQUE : 09/02/99 : ajout préconditionneur &lt;> MATRIK
C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau cgilu0
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 : 22/03/00 : ajout des préconditionneurs ILUT
C HISTORIQUE : 10/02/06 : Nettoyage
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
*     .. Work Vectors for CG
      POINTEUR IR.IZA,IP.IZA,IQ.IZA,IZ.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.NE.0) THEN
         IF (IMPR.GT.0) THEN
            WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
     $           'use Bi-CGSTAB instead !'
*         ENDIF
*         IRET=-1
*         GOTO 9999
         ENDIF
      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
         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
         ELSE
            WRITE(IOIMP,*) 'Erreur de programmation'
            GOTO 9999
         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 le gradient conjugué
C
      NBVA=NTT
      SEGINI IR,IP,IQ,IZ
C
C  Appel du gradient conjugué
C
      CALL CG2(KMORS,KISA,KS2B,INCX,
     $     KPREC,INVDIA,ILUM,ILUI,
     $     LRES,LNMV,
     $     IR,IP,IQ,IZ,
     $     ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
* Gestion du CTRL-C
      if (ierr.NE.0) return
C
C     Désactivation
C
      SEGSUP IR,IP,IQ,IZ
      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
*
*     Normal termination
*
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in prcg2.eso'
      RETURN
*
*     End of PRCG2.
*
      END