gmres.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\gmres.f
00002 !
00160                      SUBROUTINE GMRES
00161 !                    ****************
00162 !
00163      & (X,A,B,MESH,R0,V,AV,CFG,INFOGR,AUX)
00164 !
00165 !***********************************************************************
00166 ! BIEF   V7P0                                   12/06/2014
00167 !***********************************************************************
00168 !
00169 !
00170 !
00171 !
00172 !
00173 !
00174 !
00175 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00176 !| A              |-->| MATRIX OF THE SYSTEM
00177 !| AUX            |-->| MATRIX FOR PRECONDITIONING.
00178 !| AV             |<->| WORK ARRAY
00179 !| B              |-->| RIGHT-HAND SIDE OF THE SYSTEM
00180 !| CFG            |-->| STRUCTURE OF SOLVER CONFIGURATION
00181 !| INFOGR         |-->| IF YES, PRINT A LOG.
00182 !| MESH           |-->| MESH STRUCTURE.
00183 !| R              |<->| RESIDUAL
00184 !| V              |<->| WORK ARRAY
00185 !| X              |<->| INITIAL VALUE, THEN SOLUTION
00186 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00187 !
00188       USE BIEF, EX_GMRES => GMRES
00189 !
00190       IMPLICIT NONE
00191       INTEGER LNG,LU
00192       COMMON/INFO/LNG,LU
00193 !
00194 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00195 !
00196       TYPE(SLVCFG), INTENT(INOUT)    :: CFG
00197       TYPE(BIEF_OBJ), INTENT(INOUT)  :: B
00198       TYPE(BIEF_OBJ), INTENT(INOUT)  :: X,V,AV,R0
00199       TYPE(BIEF_MESH), INTENT(INOUT) :: MESH
00200       TYPE(BIEF_OBJ), INTENT(IN)     :: A
00201       TYPE(BIEF_OBJ), INTENT(INOUT)  :: AUX
00202       LOGICAL, INTENT(IN)            :: INFOGR
00203 !
00204 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00205 !
00206 !  MAXIMUM CFG%KRYLOV=20
00207 !
00208       DOUBLE PRECISION H(21,20),C(20),S(20),E(21),BID
00209 !
00210       DOUBLE PRECISION R,ZZ,NB,PREC,NR0
00211 !
00212       INTEGER I,J,K,L,M
00213 !
00214       LOGICAL RELAT,CROUT,GSEB,PRECO
00215 !
00216       INTRINSIC SQRT,ABS
00217 !
00218 !-----------------------------------------------------------------------
00219 !
00220       K = CFG%KRYLOV
00221 !
00222       CROUT=.FALSE.
00223       IF(MOD(CFG%PRECON,7).EQ.0) CROUT=.TRUE.
00224       GSEB=.FALSE.
00225       IF(MOD(CFG%PRECON,11).EQ.0) GSEB=.TRUE.
00226       PRECO=.FALSE.
00227       IF(CROUT.OR.GSEB.OR.MOD(CFG%PRECON,13).EQ.0) PRECO=.TRUE.
00228 !
00229       IF(PRECO) THEN
00230 !                  -1
00231 !       COMPUTES  L   B
00232         CALL GODOWN(B, AUX,B,'D',MESH,.FALSE.)
00233       ENDIF
00234 !
00235 ! INITIALISES
00236 !
00237       NB = P_DOTS(B,B,MESH)
00238       NB = SQRT(NB)
00239 !
00240       RELAT = .TRUE.
00241       IF(NB.LT.1.D0) THEN
00242         NB = 1.D0
00243         RELAT = .FALSE.
00244       ENDIF
00245 !
00246       M = 0
00247 !
00248 ! INITIALISES THE RESIDUAL R : A X0 - B
00249 !
00250       CALL MATRBL( 'X=AY    ',R0,A,X,BID,MESH)
00251 !
00252       IF(PRECO) THEN
00253         CALL GODOWN(R0, AUX,R0,'D',MESH,.FALSE.)
00254         CALL PUOG(X,AUX,X, 'D',MESH,.FALSE.)
00255       ENDIF
00256 !
00257       CALL OS( 'X=X-Y   ' , R0 , B , B , BID )
00258 !
00259 ! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
00260 !
00261       NR0=P_DOTS(R0,R0,MESH)
00262       NR0=SQRT(NR0)
00263       PREC = NR0/NB
00264 !
00265       IF(PREC.LE.CFG%EPS) GO TO 3000
00266 !
00267 !-----------------------------------------------------------------------
00268 !                   ITERATIONS LOOP
00269 !-----------------------------------------------------------------------
00270 !
00271 20    CONTINUE
00272 !
00273       M = M+1
00274 !
00275 ! COMPUTES THE VECTOR V1 = - R0 / NORM(R0)
00276 ! (- SIGN BECAUSE R = AX - B INSTEAD OF B - AX)
00277 !
00278       CALL OS('X=CY    ', V%ADR(1)%P , R0 , R0 , -1.D0/NR0 )
00279 !
00280 !-----------------------------------------------------------------------
00281 !         COMPUTES THE ORTHONORMAL BASE AND MATRIX H
00282 !-----------------------------------------------------------------------
00283 !
00284 !     K-1 1ST COLUMNS
00285 !
00286       DO J=1,K-1
00287 !
00288         IF(PRECO) THEN
00289           CALL GOUP(B , AUX , V%ADR(J)%P ,'D',MESH,.TRUE.)
00290           CALL MATRBL( 'X=AY    ',AV%ADR(J)%P,A,B,BID,MESH)
00291           CALL GODOWN(AV%ADR(J)%P,AUX,AV%ADR(J)%P,'D',MESH,.FALSE.)
00292         ELSE
00293           CALL MATRBL( 'X=AY    ',AV%ADR(J)%P,A,V%ADR(J)%P,BID,MESH)
00294         ENDIF
00295 !
00296         CALL OS('X=Y     ', V%ADR(J+1)%P ,AV%ADR(J)%P , X , BID )
00297 !
00298         DO I = 1,J
00299 !
00300           H(I,J) = P_DOTS( V%ADR(J+1)%P , V%ADR(I)%P , MESH )
00301           CALL OS('X=X+CY  ',V%ADR(J+1)%P,V%ADR(I)%P,V%ADR(I)%P,-H(I,J))
00302 !
00303         ENDDO
00304 !
00305         H(J+1,J)=P_DOTS(V%ADR(J+1)%P,V%ADR(J+1)%P,MESH)
00306         H(J+1,J) = SQRT( H(J+1,J) )
00307 !
00308         IF(H(J+1,J).LT.1.D-20) THEN
00309           IF(LNG.EQ.1) THEN
00310             WRITE(LU,*) 'GMRES : ECHEC DE L''ALGORITHME, MATRICE'
00311             WRITE(LU,*) 'DESORMAIS SUPPOSEE DIAGONALE'
00312           ENDIF
00313           IF(LNG.EQ.2) THEN
00314             WRITE(LU,*) 'GMRES: ALGORITHM FAILS, MATRIX ASSUMED'
00315             WRITE(LU,*) 'DIAGONAL'
00316           ENDIF
00317 !         SIMPLE SOLUTION IF A DIAGONAL...
00318 !         HERE NOT DONE FOR BLOCKS...
00319           IF(A%TYPE.EQ.3) THEN
00320             CALL OS('X=Y/Z   ',X=X%ADR(1)%P,Y=B%ADR(1)%P,Z=A%D)
00321           ELSE
00322             WRITE(LU,*) 'CASE NOT IMPLEMENTED IN GMRES'
00323             CALL PLANTE(1)
00324             STOP
00325           ENDIF
00326           CALL MATRBL('X=AY    ',R0,A,X,BID,MESH)
00327           CALL OS('X=X-Y   ', R0 , B , B , BID )
00328           NR0=P_DOTS(R0,R0,MESH)
00329           NR0=SQRT(NR0)
00330           PREC = NR0/NB
00331           IF(PREC.LE.CFG%EPS) THEN
00332             GO TO 3000
00333           ELSE
00334             WRITE(LU,*) 'PREC=',PREC,' CFG%EPS=',CFG%EPS
00335             WRITE(LU,*) 'GMRES: MATRIX NOT DIAGONAL,'
00336             WRITE(LU,*) '       ALGORITHM FAILS.'
00337             CALL PLANTE(1)
00338             STOP
00339           ENDIF
00340         ENDIF
00341         CALL OS('X=CX    ',V%ADR(J+1)%P, B, B, 1.D0/H(J+1,J))
00342 !
00343       ENDDO ! J
00344 !
00345 ! K-TH COLUMN (VECTOR V(K+1) IS NOT COMPLETELY BUILT)
00346 !
00347         IF(PRECO) THEN
00348           CALL GOUP(B , AUX , V%ADR(K)%P , 'D' ,MESH,.TRUE.)
00349           CALL MATRBL( 'X=AY    ',AV%ADR(K)%P,A,B,BID,MESH)
00350           CALL GODOWN(AV%ADR(K)%P,AUX,AV%ADR(K)%P,'D',MESH,.FALSE.)
00351         ELSE
00352           CALL MATRBL( 'X=AY    ',AV%ADR(K)%P,A,V%ADR(K)%P,BID,MESH)
00353         ENDIF
00354 !
00355         H(K+1,K) = P_DOTS( AV%ADR(K)%P , AV%ADR(K)%P , MESH )
00356 !
00357         DO I = 1,K
00358 !
00359           H(I,K) = P_DOTS( AV%ADR(K)%P , V%ADR(I)%P , MESH )
00360           H(K+1,K) = H(K+1,K) - H(I,K)**2
00361 !
00362         ENDDO ! I
00363 !       IN THEORY H(K+1,K) IS POSITIVE
00364 !       TO MACHINE ACCURACY
00365         H(K+1,K) = SQRT( ABS(H(K+1,K)) )
00366 !
00367 !-----------------------------------------------------------------------
00368 ! BUILDS GIVENS' ROTATIONS AND APPLIES TO H AND E
00369 !-----------------------------------------------------------------------
00370 !
00371 !     OTHER COMPONENTS FOR E ARE 0
00372       E(1) = NR0
00373 !
00374 !  ROTATIONS IN RANGE [1 - K]
00375 !
00376       DO I = 1 , K
00377 !
00378 !     MODIFIES COLUMN I OF H BY THE PREVIOUS ROTATIONS
00379       IF(I.GE.2) THEN
00380         DO J = 1,I-1
00381           ZZ       =  C(J) * H(J,I) + S(J) * H(J+1,I)
00382           H(J+1,I) = -S(J) * H(J,I) + C(J) * H(J+1,I)
00383           H(J,I) = ZZ
00384         ENDDO ! J
00385       ENDIF
00386 !     MODIFIES COLUMN I OF H BY ROTATION I
00387       R = SQRT( H(I,I)**2 + H(I+1,I)**2 )
00388       IF(ABS(R).LT.1.D-6) THEN
00389         IF(INFOGR) THEN
00390           IF (LNG.EQ.1) WRITE(LU,91) R
00391           IF (LNG.EQ.2) WRITE(LU,92) R
00392         ENDIF
00393         GO TO 3000
00394       ENDIF
00395       C(I) =  H(I,I)   / R
00396       S(I) =  H(I+1,I) / R
00397       H(I,I) = R
00398 !     H(I+1,I) = 0.D0    (WILL NOT BE USED AGAIN)
00399 !     MODIFIES VECTOR E
00400       E(I+1) = -S(I) * E(I)
00401       E(I  ) =  C(I) * E(I)
00402 !
00403       ENDDO ! I
00404 !
00405 !-----------------------------------------------------------------------
00406 ! SOLVES SYSTEM H*Y = E     (H UPPER TRIANGULAR OF DIMENSION K)
00407 !                            Y SAME AS E
00408 !
00409 ! H(I,I) <> 0 HAS ALREADY BEEN CHECKED ON R
00410 !-----------------------------------------------------------------------
00411 !
00412       E(K) = E(K) / H(K,K)
00413       DO J = K-1,1,-1
00414       DO L = J+1,K
00415         E(J) = E(J) - H(J,L) * E(L)
00416       ENDDO ! L
00417       E(J) = E(J) / H(J,J)
00418       ENDDO ! J
00419 !
00420 !-----------------------------------------------------------------------
00421 ! BUILDS THE SOLUTION FOR STEP M : X(M+1) = X(M) + VK * Y(K)
00422 !-----------------------------------------------------------------------
00423 !
00424       DO L = 1,K
00425 !
00426 !  COMPUTES THE NEW SOLUTION X
00427 !
00428         CALL OS('X=X+CY  ', X , V%ADR(L)%P , X , E(L) )
00429 !
00430 !  COMPUTES THE RESIDUAL : RM+1 = RM + A * VK * ZK
00431 !
00432         CALL OS('X=X+CY  ', R0 , AV%ADR(L)%P , R0 , E(L) )
00433 !
00434       ENDDO ! L
00435 !
00436 ! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
00437 ! THE RESIDUAL NORM IS GIVEN BY ABS(E(K+1))
00438 !
00439 !     NR0  =NORM OF R0
00440       NR0  = ABS(E(K+1))
00441       PREC = NR0/NB
00442 !
00443       IF (PREC.GE.CFG%EPS.AND.M.LT.CFG%NITMAX)  GOTO 20
00444 !
00445 3000  CONTINUE
00446 !
00447 !-----------------------------------------------------------------------
00448 !
00449       IF(PRECO) THEN
00450         CALL GOUP( X,AUX,X,'D',MESH,.FALSE.)
00451       ENDIF
00452 !
00453 !-----------------------------------------------------------------------
00454 !
00455 !     IF(INFOGR) THEN
00456         IF(M.LT.CFG%NITMAX) THEN
00457           IF(INFOGR) THEN
00458           IF(RELAT) THEN
00459             IF (LNG.EQ.1) WRITE(LU,101) M,PREC
00460             IF (LNG.EQ.2) WRITE(LU,102) M,PREC
00461           ELSE
00462             IF (LNG.EQ.1) WRITE(LU,201) M,PREC
00463             IF (LNG.EQ.2) WRITE(LU,202) M,PREC
00464           ENDIF
00465           ENDIF
00466         ELSE
00467           IF(RELAT) THEN
00468             IF (LNG.EQ.1) WRITE(LU,103) M,PREC
00469             IF (LNG.EQ.2) WRITE(LU,104) M,PREC
00470           ELSE
00471             IF (LNG.EQ.1) WRITE(LU,203) M,PREC
00472             IF (LNG.EQ.2) WRITE(LU,204) M,PREC
00473           ENDIF
00474         ENDIF
00475 !     ENDIF
00476 !
00477       RETURN
00478 !
00479 !-----------------------------------------------------------------------
00480 !
00481  91   FORMAT(1X,'GMRES (BIEF) : ECHEC DE L''ALGORITHME  R= ',G16.7,/,
00482      &       1X,'               SI LA MATRICE EST DIAGONALE, PRENDRE',/,
00483      &       1X,'               COMME SOLVEUR LE GRADIENT CONJUGUE.')
00484  92   FORMAT(1X,'GMRES (BIEF) : ALGORITHM FAILED  R=',G16.7,/,
00485      &       1X,'               IF THE MATRIX IS DIAGONAL, CHOOSE',/,
00486      &       1X,'               THE CONJUGATE GRADIENT SOLVER.')
00487 101   FORMAT(1X,'GMRES (BIEF) : ',
00488      &                     1I8,' ITERATIONS, PRECISION RELATIVE:',G16.7)
00489 102   FORMAT(1X,'GMRES (BIEF) : ',
00490      &                     1I8,' ITERATIONS, RELATIVE PRECISION:',G16.7)
00491 201   FORMAT(1X,'GMRES (BIEF) : ',
00492      &                     1I8,' ITERATIONS, PRECISION ABSOLUE :',G16.7)
00493 202   FORMAT(1X,'GMRES (BIEF) : ',
00494      &                     1I8,' ITERATIONS, ABSOLUTE PRECISION:',G16.7)
00495 103   FORMAT(1X,'GMRES (BIEF) : MAX D'' ITERATIONS ATTEINT:',
00496      &                     1I8,' PRECISION RELATIVE:',G16.7)
00497 104   FORMAT(1X,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00498      &                     1I8,' RELATIVE PRECISION:',G16.7)
00499 203   FORMAT(1X,'GMRES (BIEF) : MAX D'' ITERATIONS ATTEINT:',
00500      &                     1I8,' PRECISION ABSOLUE :',G16.7)
00501 204   FORMAT(1X,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00502      &                     1I8,' ABSOLUTE PRECISION:',G16.7)
00503 !
00504 !-----------------------------------------------------------------------
00505 !
00506       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0