C PRCG      SOURCE    MAGN      09/07/30    21:20:10     6458           
      SUBROUTINE PRCG(KMORS,KISA,KS2B,MATRIK,MAPREC,
     $     LRES,INCX,ITER,RESID,KPREC,
     $     IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : PRCG
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                     <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 <> 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 :
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 CCOPTIO
-INC SMLREEL
      INTEGER JG
      POINTEUR LRES.MLREEL
-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.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.6) 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
         INVDIA=0
      ENDIF
C
C     Initialisation de l'historique de convergence
C
      JG=ITER
      SEGINI LRES
C
C
C     Initialisation des vecteurs de travail pour le gradient conjugué
C
      NBVA=NTT
      SEGINI IR,IP,IQ,IZ
C
      IF (KPREC.EQ.0) THEN
C
C     Le gradient conjugué non préconditionné (fortran pur)
C
         CALL CG(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        LRES.PROG,
     $        IR.A,IP.A,IQ.A,IZ.A,
     $        ITER,RESID,IMPR,IRET)
      ELSE IF (KPREC.EQ.1) THEN
C
C     Le gradient conjugué préconditionné par la diagonale (fortran pur)
C
         CALL CGD(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        INVDIA.A,
     $        LRES.PROG,
     $        IR.A,IP.A,IQ.A,IZ.A,
     $        ITER,RESID,IMPR,IRET)
      ELSE IF (KPREC.EQ.2) THEN
C
C     Le gradient conjugué préconditionné par une factorisation
C     incomplete D-ILU (fortran pur)
C
         CALL CGDILU(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        INVDIA.A,
     $        LRES.PROG,
     $        IR.A,IP.A,IQ.A,IZ.A,
     $        ITER,RESID,IMPR,IRET)
      ELSE IF (KPREC.GE.3.AND.KPREC.LE.6) THEN
C
C     Le gradient conjugué préconditionné par une factorisation
C     incomplete (fortran pur)
C
         CALL CGILU0(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        ILUM.IA,ILUM.JA,ILUI.A,
     $        LRES.PROG,
     $        IR.A,IP.A,IQ.A,IZ.A,
     $        ITER,RESID,IMPR,IRET)
      ELSE
         WRITE (IOIMP,*) 'KPREC out of range (0..6)'
         IRET=-3
         GOTO 9999
      ENDIF
C
C     Désactivation
C
      SEGSUP IR,IP,IQ
      JG=ITER
      SEGADJ LRES
      SEGDES LRES
      IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
         SEGDES INVDIA
      ELSEIF (KPREC.GE.3.AND.KPREC.LE.6) 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 prcg.eso'
      RETURN
*
*     End of PRCG.
*
      END