rescjg.f

Go to the documentation of this file.

!                 *****************
                  SUBROUTINE rescjg
!                 *****************
!
     &(x, a,b , mesh,d,ad,ag,g,r, cfg,infogr,aux)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    SOLVES THE LINEAR SYSTEM A X = B
!+                USING THE CONJUGATE RESIDUAL METHOD.
!code
!+-----------------------------------------------------------------------
!+                        PRECONDITIONING
!+-----------------------------------------------------------------------
!+    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 (LNHE)
!+        27/02/06
!+        V5P6
!+
!
!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.
!| AG             |<->| WORK ARRAY: MATRICE A MULTIPLIED BY G.
!| 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 bief, ex_rescjg => rescjg
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL, INTENT(IN) :: INFOGR
!
!     STRUCTURES
!
      TYPE(bief_obj), INTENT(INOUT) :: D,AD,G,AG,R,X,B
      TYPE(slvcfg)  , INTENT(INOUT) :: CFG
!
!     MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!     MATRIX STRUCTURE
!
      TYPE(bief_obj), INTENT(IN)    :: A
      TYPE(bief_obj), INTENT(INOUT) :: AUX
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER M
!
      DOUBLE PRECISION XL,RMRM,TESTL,GAD
      DOUBLE PRECISION AGAD,BETA,ADAD,RO,DAD
      DOUBLE PRECISION C
!
      LOGICAL RELAT,PREC,CROUT,GSEB,PREBE,PRE3D
!
!-----------------------------------------------------------------------
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
!   INITIALISES
!
      crout =.false.
      IF(7*(cfg%PRECON/7).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)
      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)
!
      IF (rmrm.LT.cfg%EPS**2*xl) GO TO 900
!
!-----------------------------------------------------------------------
! PRECONDITIONING :
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!
! COMPUTES C G0 = 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
!
      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)
!
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!
!   COMPUTES  C DPRIM = AD  (DPRIM PUT IN B)
!
        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(b, aux , ad , 'D' , mesh)
        ELSEIF(pre3d) THEN
          CALL cpstvc(r%ADR(1)%P,g%ADR(1)%P)
          CALL trid3d(aux%X%R,b%ADR(1)%P%R,ad%ADR(1)%P%R,
     &                mesh%NPOIN,bief_nbpts(11,mesh))
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES INITIAL RO :
!-----------------------------------------------------------------------
!
      dad = p_dots(d,ad,mesh)
      IF(prec) THEN
        adad = p_dots(ad,b,mesh)
      ELSE
        adad = p_dots(ad,ad,mesh)
      ENDIF
      ro = dad/adad
!
!-----------------------------------------------------------------------
!
! COMPUTES X1 = X0 - RO  * D
!
      CALL os('X=X+CY  ', x=x, y=d, c=-ro)
!
!-----------------------------------------------------------------------
!  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
!
      rmrm   = p_dots(r,r,mesh)
!
! CHECKS END:
!
      IF (rmrm.LE.xl*cfg%EPS**2) GO TO 900
!
!-----------------------------------------------------------------------
! PRECONDITIONING : SOLVES C G = R
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!       UPDATES G BY RECURRENCE (IN B: DPRIM)
        CALL os('X=X+CY  ', x=g, y=b, c=-ro)
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES AG :
!-----------------------------------------------------------------------
!
      CALL matrbl( 'X=AY    ',ag,a,g,  c,mesh)
!
!-----------------------------------------------------------------------
! COMPUTES D BY RECURRENCE:
!-----------------------------------------------------------------------
!
      IF(prec) THEN
        agad = p_dots(ag,b,mesh)
      ELSE
        agad = p_dots(ag,ad,mesh)
      ENDIF
      beta = - agad / adad
!
      CALL os('X=CX    ', x=d, c=beta)
      CALL os('X=X+Y   ', x=d, y=g)
!
!-----------------------------------------------------------------------
! COMPUTES A D :
!-----------------------------------------------------------------------
!
      CALL os('X=CX    ', x=ad, c=beta)
      CALL os('X=X+Y   ', x=ad, y=ag)
!
      IF(prec) THEN
!
!   COMPUTES  C DPRIM = AD  (DPRIM PUT IN B)
!
        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(b , aux , ad , 'D' , mesh)
        ELSEIF(pre3d) THEN
          CALL trid3d(aux%X%R,b%ADR(1)%P%R,ad%ADR(1)%P%R,
     &                mesh%NPOIN,bief_nbpts(11,mesh))
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES RO
!-----------------------------------------------------------------------
!
      gad = p_dots(g,ad,mesh)
      IF(prec) THEN
        adad = p_dots(ad,b,mesh)
      ELSE
        adad = p_dots(ad,ad,mesh)
      ENDIF
      ro = gad/adad
!
! COMPUTES X(M) = X(M-1) - RO * D
!
      CALL os('X=X+CY  ', x=x, y=d, c=-ro)
!
      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  RETURN
!
!-----------------------------------------------------------------------
!
!   FORMATS
!
102   FORMAT(1x,'RESCJG (BIEF) : ',
     &                     1i8,' ITERATIONS, RELATIVE PRECISION:',g16.7)
202   FORMAT(1x,'RESCJG (BIEF) : ',
     &                     1i8,' ITERATIONS, ABSOLUTE PRECISION:',g16.7)
104   FORMAT(1x,'RESCJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' RELATIVE PRECISION:',g16.7)
204   FORMAT(1x,'RESCJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' ABSOLUTE PRECISON:',g16.7)
!
!-----------------------------------------------------------------------
!
      END