C BCGS      SOURCE    CHAT      05/01/12    21:37:09     5004           
      SUBROUTINE BCGS(N,NNZ,ROWPTR,COLIND,VAL,
     $     B,X,
     $     LRES,
     $     R,RTLD,P,PHAT,
     $     V,S,SHAT,T,
     $     ITER,BRTOL,RESID,IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : BCGS
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) non 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) : DCOPY, DAXPY      (subroutines)
C                    DNRM2, DDOT,DSCAL (functions)
C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
C***********************************************************************
C ENTREES            : N, NNZ, 
C                      ROWPTR, COLIND, VAL, 
C                      B, BRTOL, IMPR
C ENTREES/SORTIES    : X, 
C                      R, RTLD, P, PHAT, 
C                      V, S, SHAT, T,
C                      ITER, RESID
C SORTIES            : LRES, IRET
C CODE RETOUR (IRET) : 0 si ok
C                     <0 si problème (cf. BRTOL)
C                     >0 si l'algorithme n'a pas convergé
C     N      : nombre de degrés de liberté
C     NNZ    : nombre de valeurs non nulles de la matrice Morse
C     ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
C                           et valeurs de la matrice Morse
C     B      : vecteur second membre
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     X      : 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     R,RTLD,P,PHAT: vecteurs de travail initialisés et détruits
C     V,S,SHAT,T     dans prbcgs.
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 :
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 CCOPTIO
      INTEGER IMPR
*     .. Scalar Arguments ..
      INTEGER  N,NNZ
      INTEGER  ITER,IRET
      REAL*8   BRTOL,RESID
*     ..
*     .. Array Arguments ..
*     .. Matrices stockées en Morse
      INTEGER ROWPTR(N+1)
      INTEGER COLIND(NNZ)
      REAL*8  VAL(NNZ)
*     .. Vecteurs
      REAL*8  B(N)
      REAL*8  X(N)
      REAL*8  LRES(ITER)
      REAL*8  R(N),RTLD(N),P(N),PHAT(N)
      REAL*8  V(N),S(N),SHAT(N),T(N)
*
*
*     .. Variables locales
*     .. Parameters
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D0, ONE=1.0D0)
*     ..
      INTEGER MAXIT
      REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
      REAL*8 RHOTOL,OMETOL
      REAL*8 DDOT,DNRM2
*     ..
*     .. External subroutines and functions..
      EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL
*     ..
*     .. Executable Statements ..
*
*
      IRET  = 0
      MAXIT = ITER
      TOL   = RESID
      IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgs.eso'
*
*     Set parameter tolerances.
*
      RHOTOL = BRTOL
      OMETOL = BRTOL
*
*     Set initial residual.
*
C                    r(0)=b
      CALL DCOPY(N,B,1,R,1)
      IF (DNRM2(N,X,1).NE.ZERO) THEN
C                    r(0)=b-Ax
         CALL DMOMV ('No transpose',N,-ONE, 
     $        NNZ,ROWPTR,COLIND,VAL,
     $        X,ONE,R)
      ENDIF
*
      ITER = 0

      BNRM2 = DNRM2(N,B,1)
      IF(IMPR.GT.5) THEN
         WRITE(IOIMP,*) '||B||=',BNRM2
      ENDIF

      IF (BNRM2.EQ.ZERO) BNRM2 = ONE
      RESID = DNRM2(N,R,1) / BNRM2
*      RESID = DNRM2(N,R,1)
      IF(IMPR.GT.5) THEN
         WRITE(IOIMP,*) 'Résidu initial=',RESID
      ENDIF

      IF (RESID.LE.TOL) GOTO 30
C                    r[tilde]=r(0)
      CALL DCOPY(N,R,1,RTLD,1)
*
   10 CONTINUE
*
*     Perform BiConjugate Gradient Stabilized iteration.
*
      ITER = ITER + 1
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'ITER=',ITER
      ENDIF
*
C                    rho(i-1)=r[tilde]'r(i-1) 
      RHO = DDOT(N,RTLD,1,R,1)
C                    if rho(i-1)=0 method fails
      IF (ABS(RHO).LT.RHOTOL) GOTO 25
*
*        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 DAXPY(N,-OMEGA,V,1,P,1)
         CALL DSCAL(N,BETA,P,1)
         CALL DAXPY(N,ONE,R,1,P,1)
      ELSE
C                    p(1)=r(0)
         CALL DCOPY(N,R,1,P,1)
      ENDIF
*
*        Compute direction adjusting vector PHAT and scalar ALPHA.
*
C        Preconditionner solve : Mp[hat]=p(i)
      CALL DCOPY(N,P,1,PHAT,1)
C                    v(i)=Ap[hat]
      CALL DMOMV ('No transpose',N,ONE, 
     $     NNZ,ROWPTR,COLIND,VAL,
     $     PHAT,ZERO,V)
C                    alpha(i)=rho(i-1)/(r[tilde]'v(i))
      ALPHA = RHO / DDOT(N,RTLD,1,V,1)
*
*        Early check for tolerance.
*
C                    s=r(i-1)-alpha(i)v(i)
      CALL DAXPY(N,-ALPHA,V,1,R,1)
      CALL DCOPY(N,R,1,S,1)
      RESID = DNRM2(N,S,1) / BNRM2
      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'Résidu 1 =',RESID
      ENDIF
      IF (RESID.LE.TOL) THEN
C                    x(i)=x(i-1)+alpha(i)p[hat]
         CALL DAXPY(N,ALPHA,PHAT,1,X,1)
         LRES(ITER)=RESID
         GOTO 30
      ELSE

*
*           Compute stabilizer vector SHAT and scalar OMEGA.
*
C        Preconditionner solve : Ms[hat]=s
         CALL DCOPY(N,S,1,SHAT,1)
C                    t=As[hat]
         CALL DMOMV ('No transpose',N,ONE, 
     $        NNZ,ROWPTR,COLIND,VAL,
     $        SHAT,ZERO,T)
C                    omega(i)=t's/t't
         OMEGA = DDOT(N,T,1,S,1) / 
     $           DDOT(N,T,1,T,1)
*
*           Compute new solution approximation vector X.
*
C                    x(i)=x(i-1)+alpha(i)p[hat]
C                               +omega(i)s[hat]
         CALL DAXPY(N,ALPHA,PHAT,1,X,1)
         CALL DAXPY(N,OMEGA,SHAT,1,X,1)
         IF(IMPR.GT.7) THEN
            WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
         ENDIF

*     
*           Compute residual R, check for tolerance.
*
C                    r(i)=s-omega(i)t
         CALL DAXPY(N,-OMEGA,T,1,R,1)
*
         RESID = DNRM2(N,R,1) / BNRM2 
         IF(IMPR.GT.6) THEN
            WRITE(IOIMP,*) 'Résidu 2 =',RESID
         ENDIF
         LRES(ITER)=RESID
         IF (RESID.LE.TOL) GOTO 30
*
*     Iteration fails.
*
         IF (ITER.EQ.MAXIT) THEN
            IRET = 1
            WRITE(IOIMP,*) 
     $           'bcgs.eso : Convergence to tolerance not achieved'
            WRITE(IOIMP,*) 'ITER=',ITER
            WRITE(IOIMP,*) 'TOL=',TOL
            WRITE(IOIMP,*) 'RESID=',RESID
            RETURN
         ENDIF
*     
      ENDIF
*
      IF (ABS(OMEGA).LT.OMETOL) GOTO 25
*
*     Do some more
*
      RHO1 = RHO
      GOTO 10
*
*     A breakdown in the method has occured.
*
   25 CONTINUE
      IF (ABS(RHO).LT.RHOTOL) THEN
         IRET=-1
      ELSEIF(ABS(OMEGA).LT.OMETOL) THEN
         IRET=-2
      ENDIF
      RETURN
*
*     Iteration successful
*
   30 CONTINUE
      IRET = 0
      RETURN
*
*     End of BCGS.
*
      END