equnor.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE equnor
!                   *****************
!
     &(x, a,b , mesh, d,ad,ag,g,r, cfg,infogr,aux)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SOLVES THE LINEAR SYSTEM A X = B
!+                USING METHODS OF THE TYPE CONJUGATE GRADIENT.
!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 (LNH)
!+        24/04/97
!+        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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| AD             |<->| WORK ARRAY: MATRICE A MULTIPLIED BY D.
!| AG             |<->| WORK ARRAY: A MULTIPLIED BY DESCENT GRADIENT
!| 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_equnor => equnor
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(slvcfg), INTENT(INOUT)    :: CFG
      TYPE(bief_obj), INTENT(INOUT)  :: B
      TYPE(bief_obj), INTENT(INOUT)  :: D,AD,G,AG,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,TESTL
      DOUBLE PRECISION BETA,ADAD,RO
      DOUBLE PRECISION TGMTGM,TG1TG1,C
!
      LOGICAL RELAT,PREC,CROUT,GSEB
!
!-----------------------------------------------------------------------
!
      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.
      prec=.false.
      IF(crout.OR.gseb.OR.13*(cfg%PRECON/13).EQ.cfg%PRECON) 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
        CALL downup(g, aux , r , 'D' , mesh)
!
!  C IS HERE CONSIDERED SYMMETRICAL,
!  SHOULD OTHERWISE SOLVE TC GPRIM = G
!
!        T -1
!         C   G   IS PUT IN B
!
        CALL downup(b , aux , g , 'T' , mesh)
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES THE DIRECTION OF INITIAL DESCENT:
!-----------------------------------------------------------------------
!
      IF(prec) THEN
        CALL matrbl( 'X=TAY   ',d,a,b,  c,mesh)
      ELSE
        CALL matrbl( 'X=TAY   ',d,a,g,  c,mesh)
      ENDIF
!
      tgmtgm = p_dots(d,d,mesh)
!
!-----------------------------------------------------------------------
! 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 AG)
          CALL downup(ag, aux , ad , 'D' , mesh)
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES INITIAL RO :
!-----------------------------------------------------------------------
!
      IF(prec) THEN
        adad = p_dots(ag,ag,mesh)
      ELSE
        adad = p_dots(ad,ad,mesh)
      ENDIF
      ro = tgmtgm/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 AG: DPRIM)
          CALL os('X=X+CY  ', x=g, y=ag, c=-ro)
!
          CALL downup(b , aux , g , 'T' , mesh)
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES D BY RECURRENCE:
!-----------------------------------------------------------------------
!
      IF(prec) THEN
!                                          T  T -1          T -1
!                               AD IS HERE  A  C  G    B IS  C   G
        CALL matrbl( 'X=TAY   ',ad,a,b,  c,mesh)
      ELSE
!                               AD IS HERE TAG
        CALL matrbl( 'X=TAY   ',ad,a,g,  c,mesh)
      ENDIF
!
      tg1tg1 = tgmtgm
      tgmtgm = p_dots(ad,ad,mesh)
      beta = tgmtgm / tg1tg1
!
      CALL os('X=CX    ', x=d, c=beta)
!
!                               AD IS HERE TAG
      CALL os('X=X+Y   ', x=d, y=ad)
!
!-----------------------------------------------------------------------
! COMPUTES A D :
!-----------------------------------------------------------------------
!
      CALL matrbl( 'X=AY    ',ad,a,d,  c,mesh)
!
      IF(prec) THEN
!
!           COMPUTES  C DPRIM = AD  (DPRIM PUT IN AG)
            CALL downup(ag , aux , ad , 'D' , mesh)
!
      ENDIF
!
!-----------------------------------------------------------------------
! COMPUTES RO
!-----------------------------------------------------------------------
!
      IF(prec) THEN
        adad = p_dots(ag,ag,mesh)
      ELSE
        adad = p_dots(ad,ad,mesh)
      ENDIF
      ro = tgmtgm/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,'EQUNOR (BIEF) : ',
     &                     1i8,' ITERATIONS, RELATIVE PRECISION:',g16.7)
202   FORMAT(1x,'EQUNOR (BIEF) : ',
     &                     1i8,' ITERATIONS, ABSOLUTE PRECISION:',g16.7)
104   FORMAT(1x,'EQUNOR (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' RELATIVE PRECISION:',g16.7)
204   FORMAT(1x,'EQUNOR (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' ABSOLUTE PRECISON:',g16.7)
!
!-----------------------------------------------------------------------
!
      END