cgsqua.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE cgsqua
!                   *****************
!
     &(x,a,b,mesh, g,g0,p,k,h,ahpk,cfg,infogr)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SOLVES THE LINEAR SYSTEM A X = B
!+                USING THE SQUARED CONJUGATE GRADIENT METHOD.
!code
!+     ALGORITHM:
!+
!+        |
!+        |   INITIALISATION
!+        |   ---------------
!+        |
!+        |    0            N
!+        |   X  VECTOR IN R , APPROXIMATION OF THE SOLUTION
!+        |
!+        |      0      0
!+        |     G  = A X  - B
!+        |
!+        |     K0 = P0 = G0
!+        |
!+        |
!+        |   ITERATIONS
!+        |   ----------
!+        |
!+        |              M     0
!+        |       M   ( K  ,  G  )
!+        |     RO =  ------------
!+        |               M    0
!+        |           (A P ,  G  )
!+        |
!+        |      M     M    M     M
!+        |     H   = K - RO * A P
!+        |
!+        |      M+1   M    M       M    M
!+        |     G   = G - RO * A ( H  + K  )
!+        |
!+        |      M+1   M      M     M    M
!+        |     X   = X   - RO * ( H  + K  )
!+        |
!+        |                 0     M+1
!+        |              ( G  , G   )
!+        |     BETA =   ------------
!+        |                 0   M
!+        |              ( G , G    )
!+        |
!+        |      M+1   M+1        M    M     M       M
!+        |     P  =  G   + 2*BETA  * H + BETA**2 * P
!+        |
!+        |      M+1   M+1        M    M
!+        |     K  =  G   +   BETA  * H
!+        |
!
!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  OTHER
!+-----------------------------------------------------------------------
!code
!+  MEANING OF IELM :
!+
!+  TYPE OF ELEMENT      NUMBER OF POINTS      CODED IN THIS SUBROUTINE
!+
!+  11 : P1 TRIANGLE            3                       YES
!+  12 : P2 TRIANGLE            6                       YES
!+  13 : P1-ISO P1 TRIANGLE     6                       YES
!+  14 : P2 TRIANGLE            7
!+  21 : Q1 QUADRILATERAL       4                       YES
!+  22 : Q2 QUADRILATERAL       8
!+  24 : Q2 QUADRILATERAL       9
!+  31 : P1 TETRAHEDRON         4                       YES
!+  32 : P2 TETRAHEDRON        10
!+  41 : MITHRIDATE PRISMS      6                       YES
!
!history  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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| AHPK           |<->| WORK STRUCTURE
!| B              |-->| RIGHT-HAND SIDE OF SYSTEM
!| CFG            |-->| CFG(1): STORAGE OF MATRIX
!|                |   | CFG(2): MATRIX VECTOR PRODUCT
!| G              |<->| GRADIENT.
!| G0             |<->| INITIAL GRADIENT
!| H              |<->| WORK STRUCTURE
!| INFOGR         |-->| IF YES, INFORMATION PRINTED
!| K              |<->| WORK STRUCTURE
!| MESH           |-->| MESH STRUCTURE
!| P              |<->| WORK STRUCTURE
!| X              |<--| INITIAL VALUE, THEN SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_cgsqua => cgsqua
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj)  , INTENT(INOUT) :: X,G,G0,P,K,H,AHPK
      TYPE(bief_obj)  , INTENT(IN   ) :: B
      TYPE(bief_obj)  , INTENT(IN)    :: A
      TYPE(slvcfg)    , INTENT(INOUT) :: CFG
      LOGICAL         , INTENT(IN)    :: INFOGR
      TYPE(bief_mesh) , INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION XL,RO,TESTL,RL,GMP1G0,BETA,GMG0,C
!
      INTEGER M
!
      LOGICAL RELAT
!
      INTRINSIC sqrt
!
!-----------------------------------------------------------------------
!
! COMPUTES THE NORM OF THE SECOND MEMBER
!
      xl = p_dots(b,b,mesh)
!
      IF(xl.LT.1.d0) THEN
        xl = 1.d0
        relat = .false.
      ELSE
        relat = .true.
      ENDIF
!
      m = 0
!
! INITIALISES G  : A X0 - B
!
      CALL matrbl( 'X=AY    ',g,a,x,c,  mesh)
!
      CALL os('X=X-Y   ', x=g, y=b)
!
! CHECKS THAT THE ACCURACY HAS NOT ALREADY BEEN REACHED
!
      rl   = p_dots(g,g,mesh)
!
      IF(rl.LT.cfg%EPS**2*xl) THEN
        testl = sqrt(rl/xl)
        IF (infogr) THEN
          IF(relat) THEN
            WRITE(lu,101) m,testl
          ELSE
            WRITE(lu,201) m,testl
          ENDIF
        ENDIF
        GOTO 1000
      ENDIF
!
! INITIALISES G0 , P , AND K
!
      CALL os('X=Y     ', x=g0, y=g)
      CALL os('X=Y     ', x=p , y=g)
      CALL os('X=Y     ', x=k , y=g)
!
      m = 1
20    CONTINUE
!
! COMPUTES AP (IN H, RECOMPUTED LATER)
!
      CALL matrbl( 'X=AY    ',h,a,p,c,  mesh)
!
! COMPUTES RO
!
      ro = p_dots(k,g0,mesh) / p_dots(h,g0,mesh)
!
! COMPUTES H+K (IN H, AP ALREADY BEING IN H)
!
      CALL os('X=CX    ', x=h, c=-ro)
      CALL os('X=X+CY  ', x=h, y=k, c=2.d0)
!
!            M+1   M       M       M
! COMPUTES  X   = X  -   RO * (H+K)     (H+K IN H)
!
      CALL os('X=X+CY  ', x=x, y=h, c=-ro)
!
! COMPUTES A(H+K)  (HERE H+K IN H, A(H+K) IN AHPK)
!
      CALL matrbl( 'X=AY    ',ahpk,a,h,c,  mesh)
!
!             M     0
! COMPUTES ( G   , G  )
!
      gmg0 = p_dots(g,g0,mesh)
!
! COMPUTES GM
!
      CALL os('X=X+CY  ', x=g, y=ahpk, c=-ro)
!
      rl   = p_dots(g,g,mesh)
      IF (rl.GT.cfg%EPS**2*xl) THEN
        IF (m.GE.cfg%NITMAX) THEN
            testl=sqrt(rl/xl)
            IF(relat) THEN
              WRITE(lu,103) m,testl
            ELSE
              WRITE(lu,203) m,testl
            ENDIF
          GOTO 1000
        ELSE
          m = m + 1
        ENDIF
      ELSE
        IF(infogr) THEN
          testl=sqrt(rl/xl)
          IF(relat) THEN
            WRITE(lu,101) m,testl
          ELSE
            WRITE(lu,201) m,testl
          ENDIF
        ENDIF
        GOTO 1000
      ENDIF
!
!             M+1   0
! COMPUTES ( G   , G  )
!
      gmp1g0 = p_dots(g,g0,mesh)
!
! COMPUTES BETA
!
      beta = gmp1g0 / gmg0
!
!           M
! COMPUTES H   (H+K IN H)
!
      CALL os('X=X-Y   ', x=h, y=k)
!
!           M+1
! COMPUTES P
!
      CALL os('X=CX    ', x=p, c=beta**2 )
      CALL os('X=X+Y   ', x=p, y=g)
      CALL os('X=X+CY  ', x=p, y=h, c=2*beta)
!
!           M+1
! COMPUTES K
!
      CALL os('X=Y     ', x=k, y=g)
      CALL os('X=X+CY  ', x=k, y=h, c=beta)
!
      GOTO 20
!
1000  RETURN
!
!-----------------------------------------------------------------------
!
!   FORMATS
!
101   FORMAT(1x,'CGSQUA (BIEF) : ',1i8,' ITERATIONS',
     &          ' RELATIVE PRECISION: ',g16.7)
201   FORMAT(1x,'CGSQUA (BIEF) : ',1i8,' ITERATIONS',
     &          ' ABSOLUTE PRECISION: ',g16.7)
103   FORMAT(1x,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1i8,
     &          ' RELATIVE PRECISION: ',g16.7)
203   FORMAT(1x,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1i8,
     &          ' ABSOLUTE PRECISION:',g16.7)
!
!-----------------------------------------------------------------------
!
      END