cgstab.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE cgstab
!                   *****************
!
     &(x, a,b , mesh, p,q,r,s,t,v, cfg,infogr,aux)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    SOLVES THE LINEAR SYSTEM A X = B
!+        USING THE SQUARED CONJUGATE GRADIENT METHOD STABILLISED.
!code
!+-----------------------------------------------------------------------
!+                        PRECONDITIONING
!+                        (ALSO SEE SOLV01)
!+-----------------------------------------------------------------------
!+    PRECON VALUE     I                  MEANING
!+-----------------------------------------------------------------------
!+        0            I  NO PRECONDITIONING
!+        2            I  DIAGONAL PRECONDITIONING WITH THE MATRIX
!+                     I  DIAGONAL
!+        3            I  DIAGONAL PRECONDITIONING WITH THE CONDENSED
!+                     I  MATRIX
!+        5            I  DIAGONAL PRECONDITIONING BUT ACCEPTS 0 OR
!+                     I  NEGATIVE VALUES ON THE DIAGONAL
!+        7            I  CROUT
!+-----------------------------------------------------------------------
!
!history  R RATKE (HANNOVER); A MALCHEREK (HANNOVER); J-M HERVOUET (LNH)    ; F  LEPEINTRE (LNH)
!+        24/04/97
!+        V5P1
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        10/06/2015
!+        V7P1
!+   CALL PARMOY removed, and stop if Crout preconditionning is asked
!+   in parallel.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| AUX            |-->| PRECONDITIONING MATRIX
!| B              |-->| RIGHT-HAND SIDE OF SYSTEM
!| CFG            |-->| CFG(1): STORAGE OF MATRIX
!|                |   | CFG(2): MATRIX VECTOR PRODUCT
!| INFOGR         |-->| IF YES, INFORMATION PRINTED
!| MESH           |-->| MESH STRUCTURE
!| P              |<->| WORK STRUCTURE
!| Q              |<->| WORK STRUCTURE
!| R              |<->| WORK STRUCTURE
!| S              |<->| WORK STRUCTURE
!| T              |<->| WORK STRUCTURE
!| V              |<->| WORK STRUCTURE
!| X              |<--| VALEUR INITIALE, PUIS SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_cgstab => cgstab
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj)  , INTENT(INOUT) :: X,P,Q,R,S,T,V
      TYPE(bief_obj)  , INTENT(IN)    :: AUX,A,B
      TYPE(bief_mesh) , INTENT(INOUT) :: MESH
      TYPE(slvcfg)    , INTENT(INOUT) :: CFG
      LOGICAL         , INTENT(IN)    :: INFOGR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION ALFA,ALFA1,BETA,BETA1,OMEG,OMEG1,OMEG2
      DOUBLE PRECISION XL,TESTL,RMRM,C
!
      INTEGER M
!
      LOGICAL RELAT,CROUT
!
      DOUBLE PRECISION, PARAMETER :: RMIN = 1.d-15
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
!   INITIALISES
      crout=.false.
      IF(7*(cfg%PRECON/7).EQ.cfg%PRECON) crout=.true.
!
!-----------------------------------------------------------------------
!   INITIALISES
!-----------------------------------------------------------------------
!
      m   = 0
!
!  NORM OF THE SECOND MEMBER TO COMPUTE THE RELATIVE ACCURACY:
!
      xl = p_dots(b,b,mesh)
      IF (xl.LT.1.d0) THEN
        xl = 1.d0
        relat = .false.
      ELSE
        relat = .true.
      ENDIF
!
! IF THE SECOND MEMBER IS 0, X=0 AND IT'S THE END
!
      IF(sqrt(xl).LT.rmin) THEN
        rmrm = 0.d0
        CALL os( 'X=0     ' , x=x )
        GOTO 900
      ENDIF
!
! COMPUTES THE INITIAL RESIDUAL AND EXITS IF RESIDUAL IS SMALL:
!
      CALL matrbl( 'X=AY    ',v,a,x,c,  mesh)
!
      CALL os( 'X=Y-Z   ' , x=r , y=b , z=v )
      rmrm   = p_dots(r,r,mesh)
!
      IF (rmrm.LT.cfg%EPS**2*xl) GO TO 900
!
!-----------------------------------------------------------------------
! PRECONDITIONING :
!-----------------------------------------------------------------------
!
      IF(crout) THEN
!       COMPUTES C R  = B
        IF(ncsize.GT.1) THEN
          WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
          CALL plante(1)
          stop
        ENDIF
        CALL downup(r, aux , b , 'D' , mesh)
      ELSE
        CALL os('X=Y     ', x=r, y=b)
      ENDIF
!
!-----------------------------------------------------------------------
! RESUMES INITIALISATIONS
!-----------------------------------------------------------------------
!
      IF(crout) THEN
        IF(ncsize.GT.1) THEN
          WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
          CALL plante(1)
          stop
        ENDIF
        CALL downup(v, aux , v , 'D' , mesh)
      ENDIF
!
      CALL os('X=X-Y   ', x=r, y=v)
      CALL os('X=Y     ', x=p, y=r)
      CALL os('X=0     ', x=v )
      CALL os('X=0     ', x=q )
!
      alfa  = 1.d0
      beta  = 1.d0
      omeg1 = 1.d0
!
!-----------------------------------------------------------------------
!  ITERATIONS:
!-----------------------------------------------------------------------
!
2     CONTINUE
      m  = m  + 1
!
      beta1 = p_dots(r,p,mesh)
      omeg2 = omeg1*beta1/beta
      omeg  = omeg2/alfa
      beta  = beta1
!
      CALL os('X=Y+CZ  ', x=q, y=r, z=q, c= omeg )
      CALL os('X=X+CY  ', x=q, y=v, z=v, c=-omeg2)
!
      CALL matrbl( 'X=AY    ',v,a,q,c,  mesh)
!
      IF(crout) THEN
        IF(ncsize.GT.1) THEN
          WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
          CALL plante(1)
          stop
        ENDIF
        CALL downup(v, aux , v , 'D' , mesh)
      ENDIF
!
      omeg1 = p_dots(p,v,mesh)
      omeg1 = beta1/omeg1
!
      CALL os('X=Y+CZ  ', x=s, y=r, z=v, c=-omeg1)
!
      CALL matrbl( 'X=AY    ',t,a,s,c,  mesh)
!
      IF(crout) THEN
        IF(ncsize.GT.1) THEN
          WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
          CALL plante(1)
          stop
        ENDIF
        CALL downup(t, aux , t , 'D' , mesh)
      ENDIF
!
      alfa  = p_dots(t,s,mesh)
      alfa1 = p_dots(t,t,mesh)
      alfa  = alfa/alfa1
!
      CALL os('X=X+CY  ', x=x, y=q, c=omeg1)
      CALL os('X=X+CY  ', x=x, y=s, c=alfa )
!
      CALL os('X=Y+CZ  ', x=r, y=s, z=t, c=-alfa )
!
      rmrm   = p_dots(r,r,mesh)
!
! TESTS CONVERGENCE :
!
      IF (rmrm.LE.xl*cfg%EPS**2) GO TO 900
!
      IF(m.LT.cfg%NITMAX) GO TO 2
!
!-----------------------------------------------------------------------
!
!     IF(INFOGR) THEN
        testl = sqrt( rmrm / xl )
        IF(relat) THEN
          WRITE(lu,104) m,testl
        ELSE
          WRITE(lu,204) m,testl
        ENDIF
!     ENDIF
      GO TO 1000
!
!-----------------------------------------------------------------------
!
900   CONTINUE
!
      IF(infogr) THEN
        testl = sqrt( rmrm / xl )
        IF(relat) THEN
          WRITE(lu,103) m,testl
        ELSE
          WRITE(lu,203) m,testl
        ENDIF
      ENDIF
!
1000  RETURN
!
!-----------------------------------------------------------------------
!
!   FORMATS
!
103   FORMAT(1x,'CGSTAB (BIEF) : ',1i8,' ITERATIONS',
     &          ' RELATIVE PRECISION: ',g16.7)
203   FORMAT(1x,'CGSTAB (BIEF) : ',1i8,' ITERATIONS',
     &          ' ABSOLUTE PRECISION: ',g16.7)
104   FORMAT(1x,'CGSTAB (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1i8,
     &          ' RELATIVE PRECISION: ',g16.7)
204   FORMAT(1x,'CGSTAB (BIEF) : EXCEEDING MAXIMUM ITERATIONS',1i8,
     &          ' ABSOLUTE PRECISION:',g16.7)
!
!-----------------------------------------------------------------------
!
      END