C CGSN      SOURCE    GOUNAND   11/10/07    21:15:11     7149           
      SUBROUTINE CGSN(KMORS,KISA,B,X,
     $     KPREC,INVDIA,ILUM,ILUI,
     $     LRES,LNMV,
     $     R,RTLD,P,PHAT,
     $     Q,QHAT,UHAT,
     $     XP,BP,
     $     ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER(I-N)
C***********************************************************************
C NOM         : CGSN
C DESCRIPTION :
C     Résolution d'un système linéaire Ax=b
C     par une méthode de gradient conjugué au carré (CGS) avec la 
C     stratégie de stabilisation de Neumaier préconditionnée ou non.
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, 10/02/06
C HISTORIQUE : v1, 10/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 CGSN
      POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
      POINTEUR Q.IZA,QHAT.IZA,UHAT.IZA
      POINTEUR U.IZA,VHAT.IZA
*  Note QHAT and U, UHAT,VHAT share the same memory space
      POINTEUR XP.IZA,BP.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
      REAL*8 RHOTOL
      REAL*8 GDOT,GNRM2,RNRM2
      REAL*8 MU,MUP
*     ..
*     .. External subroutines and functions..
      EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
*     ..
*     .. Executable Statements ..
*     
*     
      IRET  = 0
      INMV  = 0
      IWARN = 0
      MAXIT = ITER
      TOL   = RESID
*     
*     Same memory space
      U=QHAT
      VHAT=UHAT
*     
      IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgsn.eso'
*     
*     Set breakdown parameter tolerances.
*     
      RHOTOL = BRTOL
      TAUTOL = BRTOL
*     
*     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(IMVEC,'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
      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)
C                    Neumaier's Strategy
C                    x=x(0)  x'=0  
C          Normalement, c'est déjà le cas
C                    b'=r(0)   mu' = ||b'||
      CALL GCOPY(R,BP)
      MUP=RNRM2
*!      WRITE(IOIMP,*) 'MUP0=',MUP
*!      WRITE(IOIMP,*) '||XP||',GNRM2(XP)
*     
 10   CONTINUE
*     
*     Perform Preconditionned Conjugate Gradient Squared 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,IDDOT)
*      IF ((ITER.LE.10).OR.((ITER/100)*100.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))
         BETA = RHO / RHO1 
C                    u(i)=r(i-1)+beta(i-1)*q(i-1)
         CALL GCOPY(R,U)
         CALL GAXPY(BETA,Q,U)
C                    p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
         CALL GSCAL(BETA**2,P)
         CALL GAXPY(BETA,Q,P)
         CALL GAXPY(ONE,U,P)
      ELSE
C                    u(1)=r(0)
C                    p(1)=u(1)
         CALL GCOPY(R,U)
         CALL GCOPY(U,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                    vhat=Ap[hat]
      CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
      INMV=INMV+1
C                    alpha(i)=rho(i-1)/(r[tilde]'vhat)
      TAU   = GDOT(RTLD,VHAT,IDDOT)
      IF (ABS(TAU).LT.TAUTOL) THEN
         TAU=SIGN(XZPREC,TAU)
         IWARN=IWARN+1
      ENDIF
      ALPHA = RHO / TAU
C                    q(i)=u(i)-alpha(i)vhat
      CALL GCOPY(U,Q)
      CALL GAXPY(-ALPHA,VHAT,Q)
*     
*        Compute direction adjusting vector UHAT.
*        PHAT is being used as temporary storage here.
*     
C                    M uhat=u(i)+q(i)
      CALL GCOPY(Q,PHAT)
      CALL GAXPY(ONE,U,PHAT)
      CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
     $     ,UHAT,PHAT)
*     
*        Compute new solution approximation vector X.
*     
*                    x'(i)=x'(i-1)+alpha(i)uhat
      CALL GAXPY(ALPHA,UHAT,XP)
*     
*        Compute residual R and check for tolerance
*        using Neumaier's strategy   
*
C                    r(i)=b'-Ax'(i)
      CALL GCOPY(BP,R)
      CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
      INMV=INMV+1
      MU=GNRM2(R)
      RESID = MU / BNRM2
CC                    qhat= A uhat
C      CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
C      INMV=INMV+1
CC                    r(i)=r(i-1)-alpha(i) qhat
C      CALL GAXPY(-ALPHA,QHAT,R)
C      RESID = GNRM2(R) / BNRM2


      IF(IMPR.GT.6) THEN
         WRITE(IOIMP,*) 'Résidu =',RESID
      ENDIF
      LRES.PROG(ITER+1)=RESID
      LNMV.LECT(ITER+1)=INMV
      IF (RESID.LE.TOL) THEN
*                    x(i)=x(i)+x'(i)         
         CALL GAXPY(ONE,XP,X)
         GOTO 30
      ENDIF
*        Neumaier's strategy   
      IF (MU.LT.MUP) THEN
*                    x(i)=x(i)+x'(i)         
         CALL GAXPY(ONE,XP,X)
*                    x'(i)=0
         CALL GSCAL(ZERO,XP)
*                    b'=r(i)
         CALL GCOPY(R,BP)
*                    mu'=mu
         MUP = MU
*         WRITE(IOIMP,*) 'ITER=',ITER
*         WRITE(IOIMP,*) 'MUP=',MUP
      ENDIF
*     
*     Iteration fails.
*     
*      IF (ITER.EQ.MAXIT) THEN
      IF (INMV.GE.MAXIT) THEN
         IRET = 1
         IF (IMPR.GT.0) THEN
C            WRITE(IOIMP,*) 
C     $           'cgsn.eso : Convergence to tolerance not achieved'
C            WRITE(IOIMP,*) 'INMV=',INMV,
C     $           ' TOL=',TOL,
C     $           ' RESID=',RESID
            IF (IWARN.GT.0) THEN
               WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
     $              ,' times.'
            ENDIF
         ENDIF
         RETURN
      ENDIF
*     
*     Do some more
*     
      RHO1 = RHO
      GOTO 10
*     
*     Iteration successful
*     
 30   CONTINUE
      IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
         WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
      ENDIF
      IRET = 0
      RETURN
*     
*     End of CGSN.
*     
      END