gmres.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE gmres
!                   ****************
!
     & (x,a,b,mesh,r0,v,av,cfg,infogr,aux)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    SOLVES A LINEAR SYSTEM A X = B
!+                USING THE GMRES (GENERALISED MINIMUM RESIDUAL) METHOD.
!code
!+  ALGORITHM:
!+
!+        |
!+        |   INITIALISATION
!+        |   ---------------
!+        |
!+        |   R0 = B - A X0
!+        |
!+        |   V1 = R0 / IIR0II
!+        |
!+        |   K = DIMENSION OF THE KRYLOV SPACE
!+        |
!+        |
!+        |
!+        |
!+        |   ITERATIONS
!+        |   ----------
!+        |
!+        | 1) BUILDS AN ORTHONORMAL BASE (VI) OF DIMENSION K+1
!+        |
!+        |   I = 1,...,K
!+        |
!+        |       VI+1 = A * VI
!+        |
!+        |       J = 1,...,I
!+        |
!+        |       HI+1,J = ( VI+1 , VJ )   HESSENBERG MATRIX (K+1,K)
!+        |
!+        |       VI+1 <--- VI+1  -  HI+1,J * VJ
!+        |
!+        |       VI+1 = VI+1 / IIVI+1II
!+        |
!+        |
!+        | 2) MULTIPLIES H BY Q  =  RK * RK-1 * ... * R1
!+        |
!+        |  GIVENS' ROTATION MATRIX RI :
!+        |                                    -                     -
!+        |                                   I ID(J-1)               I
!+        |                                   I                       I
!+        |                                   I        CJ  SJ         I
!+        |                              RI = I       -SJ  CJ         I
!+        |                                   I                       I
!+        |                                   I              ID(K-J)  I
!+        |                                    -                     -
!+        |
!+        |                         2    2
!+        |       WITH      CJ + SJ  =  1
!+        |
!+        |                 CJ AND SJ SUCH THAT H BECOMES TRIANGULAR
!+        |
!+        |                 ID(J-1) :    IDENTITY MATRIX J-1*J-1
!+        |
!+        |
!+        | 3) SOLVES THE SYSTEM H Y = Q E
!+        |                         -      -
!+        |       E VECTOR (NR0,0,0,0,0,0,.....)
!+        |
!+        |       NR0 NORM OF THE RESIDUAL
!+        |
!+        |
!+        |
!+        | 4) COMPUTES X(M+1) = X(M)  +  V * Y
!+        |
!+        |  V : MATRIX (3*NPOIN,K) WHICH COLUMNS ARE THE VJ
!+        |
!+        | 5) CHECKS THE RESIDUAL...
!+        |
!
!warning  CROUT AND GSEBE PRECONDITIONING ARE TREATED HERE LIKE A
!+            LU PRECONDITIONING. FOR CROUT IT MEANS THAT THE DIAGONAL
!+            OF THE PRECONDITIONING MATRIX IS NOT TAKEN INTO ACCOUNT
!code
!+-----------------------------------------------------------------------
!+                        PRECONDITIONING
!+-----------------------------------------------------------------------
!+    PRECON VALUE     I                  MEANING
!+-----------------------------------------------------------------------
!+        0 OR 1       I  NO PRECONDITIONING
!+                     I
!+        2            I  DIAGONAL PRECONDITIONING USING THE MATRIX
!+                     I  DIAGONAL
!+                     I
!+        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 & C. MOULIN (LNH)
!+        26/08/1993
!+        FIRST IMPLEMENTATION
!+   First version.
!
!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)
!+        12/06/2014
!+        V7P0
!+   The case A=diagonal caused a crash, because in that case the Krylov
!+   space cannot be a basis. This case is now treated since the
!+   solution is easily found : X=A%D/B
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        23/07/2015
!+        V7P1
!+   Now when the algorithm fails the dimension of the Krylov space
!+   is reduced. This is more general than the previous solution
!+   that assumed the matrix diagonal. Various cases with diagonal
!+   matrix also tested.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM
!| AUX            |-->| MATRIX FOR PRECONDITIONING.
!| AV             |<->| WORK ARRAY
!| B              |-->| RIGHT-HAND SIDE OF THE SYSTEM
!| CFG            |-->| STRUCTURE OF SOLVER CONFIGURATION
!| INFOGR         |-->| IF YES, PRINT A LOG.
!| MESH           |-->| MESH STRUCTURE.
!| R              |<->| RESIDUAL
!| V              |<->| WORK ARRAY
!| X              |<->| INITIAL VALUE, THEN SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_gmres => gmres
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(slvcfg), INTENT(INOUT)    :: CFG
      TYPE(bief_obj), INTENT(INOUT)  :: B
      TYPE(bief_obj), INTENT(INOUT)  :: X,V,AV,R0
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
      TYPE(bief_obj), INTENT(IN)     :: A
      TYPE(bief_obj), INTENT(INOUT)  :: AUX
      LOGICAL, INTENT(IN)            :: INFOGR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     HARDCODED MAXIMUM CFG%KRYLOV=20
!
      DOUBLE PRECISION H(21,20),C(20),S(20),E(21),BID
!
      DOUBLE PRECISION R,ZZ,NB,PREC,NR0
!
      INTEGER I,J,K,L,M
!
      LOGICAL RELAT,CROUT,GSEB,PRECO
!
      INTRINSIC sqrt,abs
!
!-----------------------------------------------------------------------
!
      crout=.false.
      IF(mod(cfg%PRECON,7).EQ.0) crout=.true.
      gseb=.false.
      IF(mod(cfg%PRECON,11).EQ.0) gseb=.true.
      preco=.false.
      IF(crout.OR.gseb.OR.mod(cfg%PRECON,13).EQ.0) preco=.true.
!
      IF(preco) THEN
!                  -1
!       COMPUTES  L   B
        CALL godown(b, aux,b,'D',mesh,.false.)
      ENDIF
!
! INITIALISES
!
      nb = p_dots(b,b,mesh)
      nb = sqrt(nb)
!
      relat = .true.
      IF(nb.LT.1.d0) THEN
        nb = 1.d0
        relat = .false.
      ENDIF
!
      m = 0
!
! INITIALISES THE RESIDUAL R : A X0 - B
!
      CALL matrbl( 'X=AY    ',r0,a,x,bid,mesh)
!
      IF(preco) THEN
        CALL godown(r0, aux,r0,'D',mesh,.false.)
        CALL puog(x,aux,x, 'D',mesh,.false.)
      ENDIF
!
      CALL os('X=X-Y   ', x=r0, y=b)
!
! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
!
      nr0=p_dots(r0,r0,mesh)
      nr0=sqrt(nr0)
      prec = nr0/nb
!
      IF(prec.LE.cfg%EPS) GO TO 3000
!
!-----------------------------------------------------------------------
!                   ITERATIONS LOOP
!-----------------------------------------------------------------------
!
20    CONTINUE
!
      m = m+1
!
!     K CAN BE REDUCED DURING THE ITERATION IF THE ALGORITHM FAILS
!
      k = cfg%KRYLOV
!
! COMPUTES THE VECTOR V1 = - R0 / NORM(R0)
! (- SIGN BECAUSE R = AX - B INSTEAD OF B - AX)
!
      CALL os('X=CY    ', x=v%ADR(1)%P, y=r0, c=-1.d0/nr0)
!
!-----------------------------------------------------------------------
!         COMPUTES THE ORTHONORMAL BASE AND MATRIX H
!-----------------------------------------------------------------------
!
!     K-1 1ST COLUMNS
!
      DO j=1,k-1
!
        IF(preco) THEN
          CALL goup(b , aux , v%ADR(j)%P ,'D',mesh,.true.)
          CALL matrbl( 'X=AY    ',av%ADR(j)%P,a,b,bid,mesh)
          CALL godown(av%ADR(j)%P,aux,av%ADR(j)%P,'D',mesh,.false.)
        ELSE
          CALL matrbl( 'X=AY    ',av%ADR(j)%P,a,v%ADR(j)%P,bid,mesh)
        ENDIF
!
!       NEXT VECTOR IN THE BASIS : V%ADR(J+1)%P
!
        CALL os('X=Y     ',x=v%ADR(j+1)%P,y=av%ADR(j)%P)
!
!       ORTHOGONALISATION PROCESS
!
        DO i = 1,j
          h(i,j) = p_dots( v%ADR(j+1)%P , v%ADR(i)%P , mesh )
          CALL os('X=X+CY  ',x=v%ADR(j+1)%P,y=v%ADR(i)%P,c=-h(i,j))
        ENDDO
!
        h(j+1,j)=sqrt(p_dots(v%ADR(j+1)%P,v%ADR(j+1)%P,mesh))
!
        IF(h(j+1,j).LT.1.d-20) THEN
          IF(infogr) THEN
            WRITE(lu,*) 'GMRES: NORM OF VECTOR ',j+1
            WRITE(lu,*) 'IN THE BASIS EQUALS ',h(j+1,j)
            WRITE(lu,*) 'SIZE OF KRYLOV SPACE REDUCED TO ',j
          ENDIF
!         REDUCTION OF KRYLOV SPACE
          k=j
          GO TO 2000
        ENDIF
        CALL os('X=CX    ', x=v%ADR(j+1)%P, c=1.d0/h(j+1,j))
!
      ENDDO ! J
!
! K-TH COLUMN (VECTOR V(K+1) IS NOT COMPLETELY BUILT)
!
2000  CONTINUE
!
      IF(preco) THEN
        CALL goup(b , aux , v%ADR(k)%P , 'D' ,mesh,.true.)
        CALL matrbl( 'X=AY    ',av%ADR(k)%P,a,b,bid,mesh)
        CALL godown(av%ADR(k)%P,aux,av%ADR(k)%P,'D',mesh,.false.)
      ELSE
        CALL matrbl( 'X=AY    ',av%ADR(k)%P,a,v%ADR(k)%P,bid,mesh)
      ENDIF
!
      h(k+1,k) = p_dots( av%ADR(k)%P , av%ADR(k)%P , mesh )
!
      DO i = 1,k
        h(i,k) = p_dots( av%ADR(k)%P , v%ADR(i)%P , mesh )
        h(k+1,k) = h(k+1,k) - h(i,k)**2
      ENDDO
!
!     IN THEORY H(K+1,K) IS POSITIVE
!     TO MACHINE ACCURACY
      h(k+1,k) = sqrt( max(h(k+1,k),0.d0) )
!
!-----------------------------------------------------------------------
! BUILDS GIVENS' ROTATIONS AND APPLIES TO H AND E
!-----------------------------------------------------------------------
!
!     OTHER COMPONENTS FOR E ARE 0
      e(1) = nr0
!
!  ROTATIONS IN RANGE [1 - K]
!
      DO i = 1 , k
!
!       MODIFIES COLUMN I OF H BY THE PREVIOUS ROTATIONS
        IF(i.GE.2) THEN
          DO j = 1,i-1
            zz       =  c(j) * h(j,i) + s(j) * h(j+1,i)
            h(j+1,i) = -s(j) * h(j,i) + c(j) * h(j+1,i)
            h(j,i)   = zz
          ENDDO ! J
        ENDIF
!       MODIFIES COLUMN I OF H BY ROTATION I
        r = sqrt( h(i,i)**2 + h(i+1,i)**2 )
        IF(abs(r).LT.1.d-6) THEN
          IF(infogr) THEN
            WRITE(lu,92) r
          ENDIF
          GO TO 3000
        ENDIF
        c(i) =  h(i,i)   / r
        s(i) =  h(i+1,i) / r
        h(i,i) = r
!       H(I+1,I) = 0.D0    (WILL NOT BE USED AGAIN)
!       MODIFIES VECTOR E
        e(i+1) = -s(i) * e(i)
        e(i  ) =  c(i) * e(i)
!
      ENDDO ! I
!
!-----------------------------------------------------------------------
! SOLVES SYSTEM H*Y = E     (H UPPER TRIANGULAR OF DIMENSION K)
!                            Y SAME AS E
!
! H(I,I) <> 0 HAS ALREADY BEEN CHECKED ON R
!-----------------------------------------------------------------------
!
      e(k) = e(k) / h(k,k)
      DO j = k-1,1,-1
        DO l = j+1,k
          e(j) = e(j) - h(j,l) * e(l)
        ENDDO
        e(j) = e(j) / h(j,j)
      ENDDO
!
!-----------------------------------------------------------------------
! BUILDS THE SOLUTION FOR STEP M : X(M+1) = X(M) + VK * Y(K)
!-----------------------------------------------------------------------
!
      DO l = 1,k
!
!  COMPUTES THE NEW SOLUTION X
!
        CALL os('X=X+CY  ', x=x , y=v%ADR(l)%P , c=e(l))
!
!  COMPUTES THE RESIDUAL : RM+1 = RM + A * VK * ZK
!
        CALL os('X=X+CY  ', x=r0, y=av%ADR(l)%P, c= e(l))
!
      ENDDO
!
! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
! THE RESIDUAL NORM IS GIVEN BY ABS(E(K+1))
!
!     NR0  =NORM OF R0
      nr0  = abs(e(k+1))
      prec = nr0/nb
!
      IF (prec.GE.cfg%EPS.AND.m.LT.cfg%NITMAX)  GOTO 20
!
3000  CONTINUE
!
!-----------------------------------------------------------------------
!
      IF(preco) THEN
        CALL goup( x,aux,x,'D',mesh,.false.)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     IF(INFOGR) THEN
        IF(m.LT.cfg%NITMAX) THEN
          IF(infogr) THEN
          IF(relat) THEN
            WRITE(lu,102) m,prec
          ELSE
            WRITE(lu,202) m,prec
          ENDIF
          ENDIF
        ELSE
          IF(relat) THEN
            WRITE(lu,104) m,prec
          ELSE
            WRITE(lu,204) m,prec
          ENDIF
        ENDIF
!     ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
 92   FORMAT(1x,'GMRES (BIEF) : ALGORITHM FAILED  R=',g16.7,/,
     &       1x,'               IF THE MATRIX IS DIAGONAL, CHOOSE',/,
     &       1x,'               THE CONJUGATE GRADIENT SOLVER.')
102   FORMAT(1x,'GMRES (BIEF) : ',
     &                     1i8,' ITERATIONS, RELATIVE PRECISION:',g16.7)
202   FORMAT(1x,'GMRES (BIEF) : ',
     &                     1i8,' ITERATIONS, ABSOLUTE PRECISION:',g16.7)
104   FORMAT(1x,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' RELATIVE PRECISION:',g16.7)
204   FORMAT(1x,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
     &                     1i8,' ABSOLUTE PRECISION:',g16.7)
!
!-----------------------------------------------------------------------
!
      END