C BCS2      SOURCE    CHAT      05/01/12    21:37:40     5004           
      SUBROUTINE BCS2(N,NNZ,ROWPTR,COLIND,VAL,
     $     B,X,
     $     LRES,
     $     R,RTLD,U,
     $     V,S,W,T,
     $     ITER,BRTOL,RESID,IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : BCS2
C DESCRIPTION :
C     Résolution d'un système linéaire Ax=b
C     par la méthode du bi-gradient conjugué stabilisé(2)
C     (Bi-CGSTAB(2)) 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@InProceedings{,
C  author =       {H.A. VAN DER VORST},
C  title =        {Parallel linear systems solvers :
C                  sparse iterative methods},
C  booktitle =    {High Performance Computing in Fluid Dynamics},
C  pages =     {173-200},
C  year =      {1996},
C  editor =    {P. WESSELING},
C  volume =    {3},
C  series =    {ERCOFTAC},
C  month =     {June},
C  publisher = {Kluwer Academic Publishers},
C}
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, U,
C                      V, S, W, 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,U: vecteurs de travail initialisés et détruits
C     V,S,W,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 :     15/06/98, correction d'un bug sur le calcul
C                            de tau
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),U(N)
      REAL*8  V(N),S(N),W(N),T(N)
*
*
*     .. Variables locales
*     .. Parameters
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D0, ONE=1.0D0)
*     ..
      INTEGER MAXIT, III
      REAL*8 TOL,ALPHA,BETA,RHO0,RHO1,BNRM2
      REAL*8 MU,NU,TAU,OMEGA1,OMEGA2
      REAL*8 RHOTOL
      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 bcs2.eso'
*
*     Set parameter tolerances.
*
      RHOTOL = 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.6) 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.6) 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)
*
      RHO0=1
      ALPHA=0
      OMEGA2=1
*
   10 CONTINUE
*
*     Perform BiConjugate Gradient Stabilized(2) iteration.
*
      ITER = ITER + 1
      IF(IMPR.GT.7) THEN
         WRITE(IOIMP,*) 'ITER=',ITER
      ENDIF
*
      RHO0=-OMEGA2*RHO0
*
*        Even BiCG step
*
      RHO1 = DDOT(N,RTLD,1,R,1)
      IF (ABS(RHO1).LT.RHOTOL) GOTO 25
      BETA = ALPHA * ( RHO1 / RHO0 )
      RHO0 = RHO1
*
C                    u=r(i)-beta*u
      CALL DSCAL(N,-BETA,U,1)
      CALL DAXPY(N,ONE,R,1,U,1)
C                    v=Au
      CALL DMOMV ('No transpose',N,ONE, 
     $     NNZ,ROWPTR,COLIND,VAL,
     $     U,ZERO,V)
C                    alpha=rho(0)/(r[tilde]'v)
      ALPHA = RHO0 / DDOT(N,RTLD,1,V,1)
C                    r=r(i)-alpha*v
      CALL DAXPY(N,-ALPHA,V,1,R,1)
C                    s=Ar
      CALL DMOMV ('No transpose',N,ONE, 
     $     NNZ,ROWPTR,COLIND,VAL,
     $     R,ZERO,S)
C                    x=x(i)+alpha*u
      CALL DAXPY(N,ALPHA,U,1,X,1)
*
*        Odd BiCG step
*
      RHO1 = DDOT(N,RTLD,1,S,1)
      IF (ABS(RHO1).LT.RHOTOL) GOTO 25
      BETA = ALPHA * ( RHO1 / RHO0 )
      RHO0 = RHO1
*
C                    v=s-beta*v
      CALL DSCAL(N,-BETA,V,1)
      CALL DAXPY(N,ONE,S,1,V,1)
C                    w=Av
      CALL DMOMV ('No transpose',N,ONE, 
     $     NNZ,ROWPTR,COLIND,VAL,
     $     V,ZERO,W)
C                    alpha=rho(0)/(r[tilde]'w)
      ALPHA = RHO0 / DDOT(N,RTLD,1,W,1)
C                    u=r-beta*u
      CALL DSCAL(N,-BETA,U,1)
      CALL DAXPY(N,ONE,R,1,U,1)
C                    r=r-alpha*v
      CALL DAXPY(N,-ALPHA,V,1,R,1)
C                    s=s-alpha*w
      CALL DAXPY(N,-ALPHA,W,1,S,1)
C                    t=As
      CALL DMOMV ('No transpose',N,ONE, 
     $     NNZ,ROWPTR,COLIND,VAL,
     $     S,ZERO,T)
*
*        GCR(2)-part
*
      OMEGA1=DDOT(N,R,1,S,1)
      MU    =DDOT(N,S,1,S,1)
      NU    =DDOT(N,S,1,T,1)
      TAU   =DDOT(N,T,1,T,1)
      OMEGA2=DDOT(N,R,1,T,1)
*
      TAU   =TAU-((NU*NU)/MU)
      OMEGA2=(OMEGA2-((NU*OMEGA1)/MU))/TAU
      OMEGA1=(OMEGA1-(NU*OMEGA2))/MU
C                    x(i+2)=x+omega(1)*r+omega(2)*s+alpha*u
      CALL DAXPY(N,OMEGA1,R,1,X,1)
      CALL DAXPY(N,OMEGA2,S,1,X,1)
      CALL DAXPY(N,ALPHA,U,1,X,1)
      IF(IMPR.GT.8) THEN
         WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
      ENDIF
C                    r(i+2)=r-omega(1)*s-omega(2)*t
      CALL DAXPY(N,-OMEGA1,S,1,R,1)
      CALL DAXPY(N,-OMEGA2,T,1,R,1)
      RESID = DNRM2(N,R,1) / BNRM2 
      IF(IMPR.GT.7) THEN
         WRITE(IOIMP,*) 'Résidu =',RESID
      ENDIF
      LRES(ITER)=RESID
      IF (RESID.LE.TOL) GOTO 30
*
*     Iteration fails.
*
      IF (ITER.EQ.MAXIT) THEN
         IRET = 1
         WRITE(IOIMP,*) 
     $        'bcs2.eso : Convergence to tolerance not achieved'
         WRITE(IOIMP,*) 'ITER=',ITER
         WRITE(IOIMP,*) 'TOL=',TOL
         WRITE(IOIMP,*) 'RESID=',RESID
         RETURN
      ENDIF
*     
*     Do some more
*
C                    u=u-omega(1)*v-omega(2)*w
      CALL DAXPY(N,-OMEGA1,V,1,U,1)
      CALL DAXPY(N,-OMEGA2,W,1,U,1)
      GOTO 10
*
*     A breakdown in the method has occured.
*
   25 CONTINUE
      IF(ABS(RHO1).LT.RHOTOL) THEN
         IRET=-2
      ENDIF
      RETURN
*
*     Iteration successful
*
   30 CONTINUE
      IRET = 0
      RETURN
*
*     End of BCS2.
*
      END