C BCG2      SOURCE    GOUNAND   11/07/21    21:15:01     7046           
      SUBROUTINE BCG2(KMORS,KISA,B,X,
     $     KPREC,INVDIA,ILUM,ILUI,
     $     LRES,LNMV,
     $     R,RTLD,P,PHAT,
     $     V,S,SHAT,T,
     $     ITER,INMV,BRTOL,RESID,ICALRS,IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
C***********************************************************************
C NOM         : BCG2
C DESCRIPTION :
C     Résolution d'un système linéaire Ax=b
C     par la méthode du bi-gradient conjugué stabilisé
C     (Bi-CGSTAB) préconditionné.
C     
C LANGAGE     : FORTRAN 77 + chouhia ESOPE (pour les E/S)
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 (BLAS 1) : GCOPY, GAXPY (subroutines)
C                    GNRM2, GDOT  (functions)
C APPELES (Calcul) : GMOMV  : produit Matrice Morse - Vecteur
C                    PSOLVE : Preconditionner solve
C***********************************************************************
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            :
C CODE RETOUR (IRET) : 0 si ok
C                     <0 si problème (non utilisé !)
C                     >0 si l'algorithme n'a pas convergé
C***********************************************************************
C VERSION    : v1, 09/02/06
C HISTORIQUE : v1, 09/02/06, création
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 .. a commenter si mise au point ok
-INC CCREEL
-INC CCOPTIO
-INC SMLREEL
      POINTEUR LRES.MLREEL
-INC SMLENTI
      POINTEUR LNMV.MLENTI
-INC SMMATRIK
      POINTEUR KMORS.PMORS
      POINTEUR KISA.IZA
      POINTEUR INVDIA.IZA
      POINTEUR ILUM.PMORS
      POINTEUR ILUI.IZA
      POINTEUR B.IZA
      POINTEUR X.IZA
*     .. Work Vectors for BiCGStab
      POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
      POINTEUR V.IZA,S.IZA,SHAT.IZA,T.IZA
*     
      INTEGER IMPR
*     .. Scalar Arguments ..
      INTEGER  N,NNZ
      INTEGER  ITER,IRET
      REAL*8   BRTOL,RESID
*     ..
*     
*     
*     .. Variables locales
*     .. Parameters
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D0, ONE=1.0D0)
*     ..
      INTEGER MAXIT,III
      REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA,TAU
      REAL*8 RHOTOL,OMETOL,TAUTOL
      REAL*8 GDOT,GNRM2
*     ..
*     .. External subroutines and functions..
      EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
*     ..
*     .. Executable Statements ..
*     
*     
      IRET  = 0
      INMV  = 0
      IWARN = 0
      IGRAN = 0
      MGRAN = 5
      MAXIT = ITER
      TOL   = RESID
      IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcg2.eso'
*     
*     Set parameter tolerances.
*     
      RHOTOL = BRTOL
      OMETOL = BRTOL
      TAUTOL = BRTOL
      RRMAX = 1.D0/(XZPREC**0.25D0)
*      RRMAX = 10.D0
*     
*     Set initial residual.
*     
C                    r(0)=b
      CALL GCOPY(B,R)
C                    r(0)=b-Ax
*      IF (GNRM2(X).NE.ZERO) THEN
C                    r(0)=b-Ax
*         CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
*      ENDIF
      CALL GMOMV('N',-ONE,KMORS,KISA,X,ONE,R)
      INMV=INMV+1
*     
      ITER = 0
      BNRM2 = GNRM2(B)
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) '||B||=',BNRM2
      ENDIF
      RNRM2 = GNRM2(R)
      IF (ICALRS.EQ.1) BNRM2=RNRM2
*     
      IF (BNRM2.LT.XPETIT) BNRM2 = ONE
      RESID = RNRM2 / BNRM2
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'Résidu initial=',RESID
      ENDIF
      RESMAX = RESID * RRMAX
      LRES.PROG(ITER+1)=RESID
      LNMV.LECT(ITER+1)=INMV
      IF (RESID.LE.TOL) GOTO 30
C                    r[tilde]=r(0)
      CALL GCOPY(R,RTLD)
* Autre choix
*      DO IN=1,RTLD.A(/1)
*         RTLD.A(IN) = ONE
*      ENDDO
*     
 10   CONTINUE
*     
*     Perform Preconditionned BiConjugate Gradient Stabilized iteration.
*     
      ITER = ITER + 1
      IF(IMPR.GT.7) THEN
         WRITE(IOIMP,*) 'ITER=',ITER
      ENDIF
*     
C                    rho(i-1)=r[tilde]'r(i-1)
      RHO = GDOT(RTLD,R)
*      IF ((ITER.LE.10).OR.((ITER/10)*10.EQ.ITER)) THEN
*         WRITE(IOIMP,*) 'ITER=',ITER,'  RHO=',RHO
*      ENDIF
C                    if rho(i-1)=0 method fails
      IF (ABS(RHO).LT.RHOTOL) THEN
         RHO=SIGN(XZPREC,RHO)
         IWARN=IWARN+1
      ENDIF
*     
*        Compute vector P.
*     
      IF (ITER.GT.1) THEN
C                    beta(i-1)=(rho(i-1)/rho(i-2))*
C                              (alpha(i-1)/omega(i-1))
         BETA = ( RHO / RHO1 ) * ( ALPHA / OMEGA )
C                    p(i)=r(i-1)+beta(i-1)*
C                         (p(i-1)-omega(i-1)*v(i-1))
         CALL GAXPY(-OMEGA,V,P)
         CALL GSCAL(BETA,P)
         CALL GAXPY(ONE,R,P)
      ELSE
C                    p(1)=r(0)
         CALL GCOPY(R,P)
      ENDIF
*     
*        Compute direction adjusting vector PHAT and scalar ALPHA.
*     
C     
C        Preconditionner solve
C                    Mp[hat]=p(i)
      CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
     $     ,PHAT,P)
      
C                    v(i)=Ap[hat]
      CALL GMOMV('N',ONE,KMORS,KISA,PHAT,ZERO,V)
      INMV=INMV+1
C                    alpha(i)=rho(i-1)/(r[tilde]'v(i))
      TAU   = GDOT(RTLD,V)
      IF (ABS(TAU).LT.TAUTOL) THEN
         TAU=SIGN(XZPREC,TAU)
         IWARN=IWARN+1
      ENDIF
      ALPHA = RHO / TAU
*     
*        Early check for tolerance.
*     
C                    s=r(i-1)-alpha(i)v(i)
      CALL GAXPY(-ALPHA,V,R)
      CALL GCOPY(R,S)
      RESID = GNRM2(S) / BNRM2
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'Résidu 1 =',RESID
      ENDIF
      IF (RESID.GT.RESMAX) THEN
         IGRAN=IGRAN+1
      ELSE
         IGRAN=0
      ENDIF
      IF (RESID.LE.TOL) THEN
C                    x(i)=x(i-1)+alpha(i)p[hat]
         CALL GAXPY(ALPHA,PHAT,X)
         LRES.PROG(ITER+1)=RESID
         LNMV.LECT(ITER+1)=INMV
         GOTO 30
      ELSE
*     
*           Compute stabilizer vector SHAT and scalar OMEGA.
*     
C     
C        Preconditionner solve (ILU(0))
C                    Ms[hat]=s(i)
         CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
     $        ,SHAT,S)
C                    t=As[hat]
         CALL GMOMV('N',ONE,KMORS,KISA,SHAT,ZERO,T)
         INMV=INMV+1
C     
C Il devrait y avoir un test ici ???
C     
C                    omega(i)=t's/t't
         OMEGA = GDOT(T,S) / GDOT(T,T)
*     
*           Compute new solution approximation vector X.
*     
C                    x(i)=x(i-1)+alpha(i)p[hat]
C                               +omega(i)s[hat]
         CALL GAXPY(ALPHA,PHAT,X)
         CALL GAXPY(OMEGA,SHAT,X)
*     
*           Compute residual R, check for tolerance.
*     
C                    r(i)=s-omega(i)t
         CALL GAXPY(-OMEGA,T,R)
*     
         RESID = GNRM2(R) / BNRM2 
         IF(IMPR.GT.6) THEN
            WRITE(IOIMP,*) 'Résidu 2 =',RESID
         ENDIF
         IF (RESID.GT.RESMAX) THEN
            IGRAN=IGRAN+1
         ELSE
            IGRAN=0
         ENDIF
         LRES.PROG(ITER+1)=RESID
         LNMV.LECT(ITER+1)=INMV
         IF (RESID.LE.TOL) GOTO 30
*     
*     Iteration fails.
*     
*         IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
         IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
            IRET = 1
            IF (IMPR.GT.0) THEN
C               WRITE(IOIMP,*) 
C     $              'bcg2.eso : Convergence to tolerance not achieved'
C               WRITE(IOIMP,*) 'INMV=',INMV,
C     $              ' TOL=',TOL,
C     $              ' RESID=',RESID
               IF (IWARN.GT.0) THEN
                  WRITE(IOIMP,*) 'BiCGStab nearly failed ',IWARN
     $                 ,' times.'
               ENDIF
            ENDIF
            RETURN
         ENDIF
*     
      ENDIF
*     
      IF (ABS(OMEGA).LT.OMETOL) THEN
         OMEGA=SIGN(XZPREC,OMEGA)
         IWARN=IWARN+1
      ENDIF
*     
*     Do some more
*     
      RHO1 = RHO
      GOTO 10
*     
*     Iteration successful
*     
 30   CONTINUE
      IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
         WRITE(IOIMP,*) 'BiCGStab nearly failed ',IWARN
     $        ,' times.'
      ENDIF
      IRET = 0
      RETURN
*     
*     End of BCG2.
*     
      END