C PRBCGS    SOURCE    GOUNAND   11/07/21    21:15:38     7046           
      SUBROUTINE PRBCGS(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,INCX,ITER,
     $     RESID,KPREC,
     $     BRTOL,IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : PRBCGS
C DESCRIPTION :
C     Préparation de la résolution d'un système linéaire Ax=b
C     par une méthode de bi-gradient conjugué stabilisé 
C     (Bi-CGSTAB <~> Bi-CG + GMRES(1))
C     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 (mais la méthode peut ne pas converger
C                             voir paramètre BRTOL en particulier).
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       bcgs, bcgsd, bcgsdi, bcgsi0
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) : BCGS, BCGSD, BCGSDI, BCGSI0
C***********************************************************************
C ENTREES            : MATRIK, MAPREC, KPREC, BRTOL, 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     BRTOL  : REAL*8
C             'Breakdown tolerance' pour les méthodes de type
C              BiCGSTAB. Si un certain produit scalaire de vecteurs
C              "direction" est inférieur à cette tolérance, la 
C              méthode ne peut pas continuer.
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 bcgsi0
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 : 09/04/04 : ajout des préconditionneurs ILUTP
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
*     .. Parameters
*     This is a breakdown tolerance for BiCGSTAB-type method
      REAL*8     BRTOL
*     .. Work Vectors for BCGS..
      POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
      POINTEUR IV.IZA,IS.IZA,ISHAT.IZA,IT.IZA
*STAT-INC SMSTAT
*     .. Scalar Arguments ..
      INTEGER ITER, KPREC, IMPR, IRET
      REAL*8  RESID
*     ..
*     Vars reqd for STOPTEST2
*      REAL*8   TOL, BNRM2
*     ..
*     .. External subroutines ..
*      EXTERNAL STOPTEST2
      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 !'
C         IRET=-2
C         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, ILUT2
C ,ilutp, ilutpg, ilutpg2
      ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
         SEGACT ILUM
         SEGACT ILUI
*         SEGPRT,ILUM
*         SEGPRT,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 BiCGSTAB
C
*STAT      CALL INMSTA(MSTAT,0)
      NBVA=NTT
      SEGINI IR,IRTLD,IP,IPHAT
      SEGINI IV,IS,ISHAT,IT
*STAT      CALL PRMSTA('    PRBCGS : objets temporaires',MSTAT,1)
*STAT      CALL SUMSTA(MSTAT,0)
C
      IF (KPREC.EQ.0) THEN
C
C     Le BiCGSTAB non préconditionné (fortran pur)
C
         CALL BCGS(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        LRES.PROG,
     $        IR.A,IRTLD.A,IP.A,IPHAT.A,
     $        IV.A,IS.A,ISHAT.A,IT.A,
     $        ITER,BRTOL,RESID,IMPR,IRET)
      ELSE IF (KPREC.EQ.1) THEN
C
C     Le gradient conjugué préconditionné par la diagonale (fortran pur)
C
         CALL BCGSD(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        INVDIA.A,
     $        LRES.PROG,
     $        IR.A,IRTLD.A,IP.A,IPHAT.A,
     $        IV.A,IS.A,ISHAT.A,IT.A,
     $        ITER,BRTOL,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 BCGSDI(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        INVDIA.A,
     $        LRES.PROG,
     $        IR.A,IRTLD.A,IP.A,IPHAT.A,
     $        IV.A,IS.A,ISHAT.A,IT.A,
     $        ITER,BRTOL,RESID,IMPR,IRET)
      ELSE IF (KPREC.GE.3.AND.KPREC.LE.9) THEN
C
C     Le gradient conjugué préconditionné par une factorisation
C     incomplete (fortran pur)
C
         CALL BCGSI0(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
     $        KS2B.A,INCX.A,
     $        ILUM.IA,ILUM.JA,ILUI.A,
     $        LRES.PROG,
     $        IR.A,IRTLD.A,IP.A,IPHAT.A,
     $        IV.A,IS.A,ISHAT.A,IT.A,
     $        ITER,BRTOL,RESID,IMPR,IRET)
      ELSE
         WRITE (IOIMP,*) 'KPREC out of range (0..9)'
         IRET=-4
         GOTO 9999
      ENDIF
C
C     Désactivation
C
      SEGSUP IR,IRTLD,IP,IPHAT
      SEGSUP IV,IS,ISHAT,IT
      JG=ITER
      SEGADJ LRES
      SEGDES LRES
      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 breakdown has occured in the BiCGSTAB method
C
      IF (IRET.LT.0) GOTO 9999
*
*     Normal termination
*
      RETURN
*
* Format handling
*
 1002 FORMAT(10(1X,1PE11.4))
*
* Error handling
*
 9999 CONTINUE
      WRITE(IOIMP,*) 'An error was detected in prbcgs.eso'
      RETURN
*
*     End of PRBCGS
*
      END