gracjg.F

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE gracjg
!                   *****************
!
     &(x, a,b , mesh, d,ad,g,r, cfg,infogr,aux)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    SOLVES THE LINEAR SYSTEM A X = B
!+                USING THE CONJUGATE GRADIENT METHOD.
!code
!+-----------------------------------------------------------------------
!+                        PRECONDITIONING
!+-----------------------------------------------------------------------
!+   CFG%PRECON VALUE  I                  MEANING
!+-----------------------------------------------------------------------
!+                     I
!+        0 OR 1       I  NO PRECONDITIONING
!+                     I
!+        2            I  DIAGONAL PRECONDITIONING USING THE MATRIX
!+                     I  DIAGONAL
!+        3            I  BLOCK-DIAGONAL PRECONDITIONING
!+                     I
!+        5            I  DIAGONAL PRECONDITIONING USING THE ABSOLUTE
!+                     I  VALUE OF THE MATRIX DIAGONAL
!+                     I
!+        7            I  CROUT EBE PRECONDITIONING
!+                     I
!+       11            I  GAUSS-SEIDEL EBE PRECONDITIONING
!+                     I
!+-----------------------------------------------------------------------
!
!history  J-M HERVOUET (LNH)
!+        06/10/08
!+        V5P9
!+
!
!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
!| AD             |<->| WORK ARRAY: MATRICE A MULTIPLIED BY D.
!| AUX            |-->| MATRIX FOR PRECONDITIONING.
!| B              |-->| RIGHT-HAND SIDE OF THE SYSTEM
!| CFG            |-->| STRUCTURE OF SOLVER CONFIGURATION
!| D              |<->| WORK ARRAY: DIRECTION OF DESCENT.
!| G              |<->| DESCENT GRADIENT.
!| INFOGR         |-->| IF YES, PRINT A LOG.
!| MESH           |-->| MESH STRUCTURE.
!| R              |<->| RESIDUAL (MAY BE IN THE SAME MEMORY SPACE AS
!|                |   | GRADIENT DEPENDING ON CONDITIONING)
!| X              |<->| INITIAL VALUE, THEN SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      USE declarations_telemac, ONLY : gracjg_cnt
      USE bief, ex_gracjg => gracjg
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(slvcfg), INTENT(INOUT)    :: CFG
      TYPE(bief_obj), INTENT(INOUT)  :: B
      TYPE(bief_obj), INTENT(INOUT)  :: D,AD,G,R,X
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
      TYPE(bief_obj), INTENT(IN)     :: A
      TYPE(bief_obj), INTENT(INOUT)  :: AUX
      LOGICAL, INTENT(IN)            :: INFOGR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER M
!
      DOUBLE PRECISION XL,RMRM,RMDM,RMGM,TESTL
      DOUBLE PRECISION BETA,RO,DAD,RM1GM1,STO1,C
!
      LOGICAL RELAT,PREC,CROUT,GSEB,PRE3D,PREBE
!
      DOUBLE PRECISION, PARAMETER :: RMIN = 1.d-15
!
!-----------------------------------------------------------------------
!
      LOGICAL           :: FINISH
!
#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
      INTEGER           :: DRV_ELMS
      DOUBLE PRECISION  :: DRV_ERR
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE  :: DRV0, DRV1
#endif
#if defined COMPAD_DCO_T1V
      INTEGER           :: T1V_NOD,I
#endif
!
!-----------------------------------------------------------------------
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
!   INITIALISES
!
      gracjg_cnt = gracjg_cnt + 1
#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
!      WRITE(LU,*) 'AD LINEAR SOLVER RESET DERIVATIVES :: ',
!     &     AD_LINSOLV_RESETDERIV
!      WRITE(LU,*) 'AD LINEAR SOLVER DERIVATIVE CONVERGENCE :: ',
!     &     AD_LINSOLV_DERIVATIVE_CONVERGENCE
!
! 1.  RESET INCOMING TANGENTS TO ZERO TO AVOID DEPENDENCY OF OVERALL
!     TANGENTS FROM INITIAL SOLUTIONS FOR FIRST CONJUGADE GRADIENT
!! COMPAD-DCO  BEGIN  JR2016
      IF (ad_linsolv_resetderiv ) THEN
#  if defined COMPAD_DCO_T1S
        CALL dco_t1s_set( x%ADR(1)%P%R, 0.d0, 1 )
#  endif
#  if defined COMPAD_DCO_T1V
        t1v_nod = dco_t1v_get_nod()
!        WRITE(LU,*) 'T1V  NOD  :: ',t1V_NOD
        DO i=1, t1v_nod
          CALL dco_t1v_set( x%ADR(1)%P%R, 0.d0, 1, i )
        END DO
#  endif
      END IF
!! COMPAD-DCO  original code
!!      CALL DCO_T1S_SET( X%ADR(1)%P%R, 0.D0, 1 )
!! COMPAD-DCO  END    JR2016
!     BY SETTING X TO A FIXED STARTING VALUE WE DISABLE THE RESTART
!     FROM THE PREVISOU CG SOLUTION, AND SETTING TANGENTS TO ZERO TOO
!
! 2.  PREPARE STORAGE OF GRADIENTS FOR TERMINATION TEST
      drv_elms = x%ADR(1)%P%DIM1
      ALLOCATE( drv0(drv_elms), drv1(drv_elms))
!
! 3.  STORE TANGENT OF X0
#  if defined COMPAD_DCO_T1S
      CALL dco_t1s_get( x%ADR(1)%P%R, drv0, 1 )
#  endif
#  if defined COMPAD_DCO_T1V
      CALL dco_t1v_get( x%ADR(1)%P%R, drv0, 1, 1 )
#  endif
!
#endif
!
!
      sto1 = 0.d0
      crout =.false.
      IF(07*(cfg%PRECON/07).EQ.cfg%PRECON) crout=.true.
      gseb=.false.
      IF(11*(cfg%PRECON/11).EQ.cfg%PRECON) gseb=.true.
      prebe=.false.
      IF(13*(cfg%PRECON/13).EQ.cfg%PRECON) prebe=.true.
      pre3d=.false.
      IF(17*(cfg%PRECON/17).EQ.cfg%PRECON) pre3d=.true.
      prec=.false.
      IF(crout.OR.gseb.OR.prebe.OR.pre3d) prec=.true.
!
!-----------------------------------------------------------------------
!   INITIALISES
!-----------------------------------------------------------------------
!
      m   = 0
!
!  NORMALISES THE SECOND MEMBER TO COMPUTE THE RELATIVE PRECISION:
!
      xl = p_dots(b,b,mesh)
!
!     OLIVIER BOITEAU'S (SINETICS) TEST : SECOND MEMBER TOO SMALL
!                                         ==> UNKNOWN = 0
!
      testl=sqrt(xl)
      IF(testl.LT.rmin) THEN
        CALL os('X=0     ',x=x)
        IF(infogr) THEN
          WRITE(lu,106) testl
        ENDIF
        GOTO 1000
      ENDIF
!
      IF(xl.LT.1.d0) THEN
        xl = 1.d0
        relat = .false.
      ELSE
        relat = .true.
      ENDIF
!
! COMPUTES THE INITIAL RESIDUAL AND POSSIBLY EXITS:
!
      CALL matrbl( 'X=AY    ',r,a,x,  c,mesh)
!
      CALL os('X=X-Y   ', x=r, y=b)
      rmrm   = p_dots(r,r,mesh)
      rmgm   = rmrm
!
      IF (rmrm.LT.cfg%EPS**2*xl) GO TO 900
!
!-----------------------------------------------------------------------
! PRECONDITIONING :
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!       COMPUTES C G0 = R0
        IF(crout.OR.gseb.OR.prebe) THEN
          IF(ncsize.GT.1) THEN
            WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
            CALL plante(1)
            stop
          ENDIF
          CALL downup( g, aux , r , 'D' , mesh )
        ELSEIF(pre3d) THEN
          CALL cpstvc(r%ADR(1)%P,g%ADR(1)%P)
          CALL trid3d(aux%X%R,g%ADR(1)%P%R,r%ADR(1)%P%R,
     &                mesh%NPOIN,bief_nbpts(11,mesh))
        ENDIF
!       COMPUTES RMGM AND STORES
        rmgm = p_dots(r,g,mesh)
        sto1 = rmgm
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES THE DIRECTION OF INITIAL DESCENT:
!-----------------------------------------------------------------------
!
      CALL os('X=Y     ', x=d, y=g)
!
!-----------------------------------------------------------------------
! COMPUTES THE INITIAL PRODUCT A D :
!-----------------------------------------------------------------------
!
      CALL matrbl( 'X=AY    ',ad,a,d,c,mesh)
!
!-----------------------------------------------------------------------
! COMPUTES INITIAL RO :
!-----------------------------------------------------------------------
!
      dad = p_dots(d,ad,mesh)
      ro = rmgm / dad
!     USES RMGM BECAUSE HERE D0=G0
!
!-----------------------------------------------------------------------
!
! COMPUTES X1 = X0 - RO  * D
!
      CALL os('X=X+CY  ',x=x,y=d,c=-ro)
!
!     ALGORITHMIC DIFFERENTIATION
!     TAKE TANGENT OF X1 IF WE ARE IN DCO_T1S MODE AND STORE IT
#if defined COMPAD_DCO_T1S
      CALL dco_t1s_get( x%ADR(1)%P%R, drv1, 1 )
#endif
#if defined COMPAD_DCO_T1V
      CALL dco_t1v_get( x%ADR(1)%P%R, drv1, 1, 1 )
#endif
!
!-----------------------------------------------------------------------
!  ITERATIONS LOOP:
!-----------------------------------------------------------------------
!
2     m  = m  + 1
!
!-----------------------------------------------------------------------
! COMPUTES THE RESIDUAL : R(M) = R(M-1) - RO(M-1) A D(M-1)
!-----------------------------------------------------------------------
!
      CALL os('X=X+CY  ',x=r,y=ad,c=-ro)
!
!  SOME VALUES WILL CHANGE IN CASE OF PRECONDITIONING
!
      rm1gm1 = rmgm
      rmrm   = p_dots(r,r,mesh)
      rmdm   = rmrm
      rmgm   = rmrm
!
! CHECKS END:
!
!      IF(RMRM.LE.XL*CFG%EPS**2) GO TO 900
      finish = rmrm.LE.xl*cfg%EPS**2
#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
!     CHECK DIFFERENCE BETWEEN TANGENTS OF LAST TWO STEPS
! 1.  SET TO ZERO TO AVOID ADDTIONAL TEST
      drv_err = 0.d0
! 2a. ABSOLUTE CHANGE TO LAST STEP
      drv_err  = sum( drv1 - drv0 )**2
! 2b. RELATIVE CHANGE TO LAST STEP
!      DRV_ERR = SUM( DRV0*DRV0 )
!      IF ( DRV_ERR .GT. 0.0 ) THEN
!        DRV_ERR  = SUM( (DRV1 - DRV0 )**2 ) / DRV_ERR
!      ENDIF
! 3.  TEST CONVERGENCE
      IF ( ad_linsolv_derivative_convergence ) THEN
        finish = finish .AND. (drv_err .LE. cfg%EPS**2 )
      ENDIF
!! COMPAD-DCO  original code
!!      FINISH = FINISH .AND. (DRV_ERR .LE. CFG%EPS**2 )
!#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
!      WRITE(LU,*) '  CG  CONVTEST',CNT,M,FINISH,RMRM,' :: ',
!     &     SQRT(DRV_ERR),(DRV_ERR.LE.CFG%EPS**2), SQRT(SUM(DRV1**2)),
!     &     SQRT(SUM(DRV0**2))
!#endif
#endif
      IF ( finish ) GO TO 900
!
!-----------------------------------------------------------------------
! PRECONDITIONING : SOLVES C G = R
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!
!       SOLVES C G = R
        IF(crout.OR.gseb.OR.prebe) THEN
          IF(ncsize.GT.1) THEN
            WRITE(lu,*) 'NO CROUT PRECONDITIONNING IN PARALLEL'
            CALL plante(1)
            stop
          ENDIF
          CALL downup( g, aux , r , 'D' , mesh )
        ELSEIF(pre3d) THEN
          CALL cpstvc(r%ADR(1)%P,g%ADR(1)%P)
          CALL trid3d(aux%X%R,g%ADR(1)%P%R,r%ADR(1)%P%R,
     &                mesh%NPOIN,bief_nbpts(11,mesh))
        ENDIF
!       COMPUTES RMGM AND RM1GM1
        rm1gm1 = sto1
        rmgm = p_dots(r,g,mesh)
        sto1=rmgm
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES D BY RECURRENCE:
!-----------------------------------------------------------------------
!
      beta = rmgm / rm1gm1
!     OPTIMISED EQUIVALENT OF THE 2 OS CALLS:
      CALL os( 'X=Y+CZ  ' , x=d , y=g , z=d , c=beta )
!
!-----------------------------------------------------------------------
! PRECONDITIONING :
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!       COMPUTES RMDM
        rmdm=p_dots(r,d,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES A D :
!-----------------------------------------------------------------------
!
      CALL matrbl( 'X=AY    ',ad,a,d,c,mesh)
!
!-----------------------------------------------------------------------
! COMPUTES RO
!-----------------------------------------------------------------------
!
      dad = p_dots(d,ad,mesh)
      ro = rmdm/dad
!
! COMPUTES X(M) = X(M-1) - RO * D
!
      CALL os('X=X+CY  ',x=x,y=d,c=-ro)
!
!-----------------------------------------------------------------------
#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
!     TAKE TANGENT OF X1 IF WE ARE IN DCO_T1S MODE
!     AND STORE TANGENT OF FORMER X
      drv0 = drv1
!     STORE TANGENT OF NEW X
#  if defined COMPAD_DCO_T1S
      CALL dco_t1s_get( x%ADR(1)%P%R, drv1, 1 )
#  endif
#  if defined COMPAD_DCO_T1V
      CALL dco_t1v_get( x%ADR(1)%P%R, drv1, 1, 1 )
#  endif
#endif
!-----------------------------------------------------------------------
!
      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,102) m,testl
        ELSE
          WRITE(lu,202) m,testl
        ENDIF
      ENDIF
!
1000  CONTINUE
!
!-----------------------------------------------------------------------
!
#if defined COMPAD_DCO_T1S || COMPAD_DCO_T1V
!      WRITE(LU,*) '  CG_END CONVTEST',
!     &  CNT,M,FINISH,RMRM,' :: ',
!     &  SQRT(DRV_ERR),(DRV_ERR.LE.CFG%EPS**2), SQRT(SUM(DRV1**2))
#endif
!
!-----------------------------------------------------------------------
!
      RETURN
!
!-----------------------------------------------------------------------
!
!   FORMATS
!
102   FORMAT(1x,'GRACJG (BIEF) : ',
     &                     1i8,' ITERATIONS, RELATIVE PRECISION:',g16.7)
202   FORMAT(1x,'GRACJG (BIEF) : ',
     &                     1i8,' ITERATIONS, ABSOLUTE PRECISION:',g16.7)
104   FORMAT(1x,'GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' RELATIVE PRECISION:',g16.7)
106   FORMAT(1x,'GRACJG (BIEF) : ',
     &        ' SOLUTION X=0 BECAUSE L2-NORM OF B VERY SMALL:',g16.7)
204   FORMAT(1x,'GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' ABSOLUTE PRECISION:',g16.7)
!
!-----------------------------------------------------------------------
!
      END