gracjg.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\gracjg.f
00002 !
00084                      SUBROUTINE GRACJG
00085 !                    *****************
00086 !
00087      &(X, A,B , MESH, D,AD,G,R, CFG,INFOGR,AUX)
00088 !
00089 !***********************************************************************
00090 ! BIEF   V6P1                                   21/08/2010
00091 !***********************************************************************
00092 !
00093 !
00094 !
00095 !
00096 !
00097 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00098 !| A              |-->| MATRIX OF THE SYSTEM
00099 !| AD             |<->| WORK ARRAY: MATRICE A MULTIPLIED BY D.
00100 !| AUX            |-->| MATRIX FOR PRECONDITIONING.
00101 !| B              |-->| RIGHT-HAND SIDE OF THE SYSTEM
00102 !| CFG            |-->| STRUCTURE OF SOLVER CONFIGURATION
00103 !| D              |<->| WORK ARRAY: DIRECTION OF DESCENT.
00104 !| G              |<->| DESCENT GRADIENT.
00105 !| INFOGR         |-->| IF YES, PRINT A LOG.
00106 !| MESH           |-->| MESH STRUCTURE.
00107 !| R              |<->| RESIDUAL (MAY BE IN THE SAME MEMORY SPACE AS
00108 !|                |   | GRADIENT DEPENDING ON CONDITIONING)
00109 !| X              |<->| INITIAL VALUE, THEN SOLUTION
00110 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00111 !
00112       USE BIEF, EX_GRACJG => GRACJG
00113 !
00114       IMPLICIT NONE
00115       INTEGER LNG,LU
00116       COMMON/INFO/LNG,LU
00117 !
00118 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00119 !
00120       TYPE(SLVCFG), INTENT(INOUT)    :: CFG
00121       TYPE(BIEF_OBJ), INTENT(INOUT)  :: B
00122       TYPE(BIEF_OBJ), INTENT(INOUT)  :: D,AD,G,R,X
00123       TYPE(BIEF_MESH), INTENT(INOUT) :: MESH
00124       TYPE(BIEF_OBJ), INTENT(IN)     :: A
00125       TYPE(BIEF_OBJ), INTENT(INOUT)  :: AUX
00126       LOGICAL, INTENT(IN)            :: INFOGR
00127 !
00128 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00129 !
00130       INTEGER M
00131 !
00132       DOUBLE PRECISION XL,RMRM,RMDM,RMGM,TESTL
00133       DOUBLE PRECISION BETA,RO,DAD,RM1GM1,STO1,C
00134 !
00135       LOGICAL RELAT,PREC,CROUT,GSEB,PRE3D,PREBE
00136 !
00137       DOUBLE PRECISION RMIN
00138       DATA RMIN/1.D-15/
00139 !
00140 !-----------------------------------------------------------------------
00141 !
00142       INTRINSIC SQRT
00143 !
00144 !-----------------------------------------------------------------------
00145 !
00146 !   INITIALISES
00147 !
00148       STO1 = 0.D0
00149       CROUT =.FALSE.
00150       IF(07*(CFG%PRECON/07).EQ.CFG%PRECON) CROUT=.TRUE.
00151       GSEB=.FALSE.
00152       IF(11*(CFG%PRECON/11).EQ.CFG%PRECON) GSEB=.TRUE.
00153       PREBE=.FALSE.
00154       IF(13*(CFG%PRECON/13).EQ.CFG%PRECON) PREBE=.TRUE.
00155       PRE3D=.FALSE.
00156       IF(17*(CFG%PRECON/17).EQ.CFG%PRECON) PRE3D=.TRUE.
00157       PREC=.FALSE.
00158       IF(CROUT.OR.GSEB.OR.PREBE.OR.PRE3D) PREC=.TRUE.
00159 !
00160 !-----------------------------------------------------------------------
00161 !   INITIALISES
00162 !-----------------------------------------------------------------------
00163 !
00164       M   = 0
00165 !
00166 !  NORMALISES THE SECOND MEMBER TO COMPUTE THE RELATIVE PRECISION:
00167 !
00168       XL = P_DOTS(B,B,MESH)
00169 !
00170 !     OLIVIER BOITEAU'S (SINETICS) TEST : SECOND MEMBER TOO SMALL
00171 !                                         ==> UNKNOWN = 0
00172 !
00173       TESTL=SQRT(XL)
00174       IF(TESTL.LT.RMIN) THEN
00175         CALL OS('X=0     ',X=X)
00176         IF(INFOGR) THEN
00177           IF(LNG.EQ.1) WRITE(LU,105) TESTL
00178           IF(LNG.EQ.2) WRITE(LU,106) TESTL
00179         ENDIF
00180         GOTO 1000
00181       ENDIF
00182 !
00183       IF(XL.LT.1.D0) THEN
00184         XL = 1.D0
00185         RELAT = .FALSE.
00186       ELSE
00187         RELAT = .TRUE.
00188       ENDIF
00189 !
00190 ! COMPUTES THE INITIAL RESIDUAL AND POSSIBLY EXITS:
00191 !
00192       CALL MATRBL( 'X=AY    ',R,A,X,  C,MESH)
00193 !
00194       CALL OS( 'X=X-Y   ' , R , B , B , C )
00195       RMRM   = P_DOTS(R,R,MESH)
00196       RMGM   = RMRM
00197 !
00198       IF (RMRM.LT.CFG%EPS**2*XL) GO TO 900
00199 !
00200 !-----------------------------------------------------------------------
00201 ! PRECONDITIONING :
00202 !-----------------------------------------------------------------------
00203 !
00204       IF(PREC) THEN
00205 !       COMPUTES C G0 = R0
00206         IF(CROUT.OR.GSEB.OR.PREBE) THEN
00207           CALL DOWNUP( G, AUX , R , 'D' , MESH )
00208           IF(NCSIZE.GT.1) CALL PARMOY(G,MESH)
00209         ELSEIF(PRE3D) THEN
00210           CALL CPSTVC(R%ADR(1)%P,G%ADR(1)%P)
00211           CALL TRID3D(AUX%X%R,G%ADR(1)%P%R,R%ADR(1)%P%R,
00212      &                MESH%NPOIN,BIEF_NBPTS(11,MESH))
00213         ENDIF
00214 !       COMPUTES RMGM AND STORES
00215         RMGM = P_DOTS(R,G,MESH)
00216         STO1 = RMGM
00217       ENDIF
00218 !
00219 !-----------------------------------------------------------------------
00220 ! COMPUTES THE DIRECTION OF INITIAL DESCENT:
00221 !-----------------------------------------------------------------------
00222 !
00223       CALL OS( 'X=Y     ' , D , G , G , C )
00224 !
00225 !-----------------------------------------------------------------------
00226 ! COMPUTES THE INITIAL PRODUCT A D :
00227 !-----------------------------------------------------------------------
00228 !
00229       CALL MATRBL( 'X=AY    ',AD,A,D,C,MESH)
00230 !
00231 !-----------------------------------------------------------------------
00232 ! COMPUTES INITIAL RO :
00233 !-----------------------------------------------------------------------
00234 !
00235       DAD = P_DOTS(D,AD,MESH)
00236       RO = RMGM / DAD
00237 !     RO = RMDM / DAD  (USES RMGM BECAUSE HERE D0=G0)
00238 !
00239 !-----------------------------------------------------------------------
00240 !
00241 ! COMPUTES X1 = X0 - RO  * D
00242 !
00243       CALL OS('X=X+CY  ',X=X,Y=D,C=-RO)
00244 !
00245 !-----------------------------------------------------------------------
00246 !  ITERATIONS LOOP:
00247 !-----------------------------------------------------------------------
00248 !
00249 2     M  = M  + 1
00250 !
00251 !-----------------------------------------------------------------------
00252 ! COMPUTES THE RESIDUAL : R(M) = R(M-1) - RO(M-1) A D(M-1)
00253 !-----------------------------------------------------------------------
00254 !
00255       CALL OS('X=X+CY  ',X=R,Y=AD,C=-RO)
00256 !
00257 !  SOME VALUES WILL CHANGE IN CASE OF PRECONDITIONING
00258 !
00259       RM1GM1 = RMGM
00260       RMRM   = P_DOTS(R,R,MESH)
00261       RMDM   = RMRM
00262       RMGM   = RMRM
00263 !
00264 ! CHECKS END:
00265 !
00266       IF(RMRM.LE.XL*CFG%EPS**2) GO TO 900
00267 !
00268 !-----------------------------------------------------------------------
00269 ! PRECONDITIONING : SOLVES C G = R
00270 !-----------------------------------------------------------------------
00271 !
00272       IF(PREC) THEN
00273 !
00274 !       SOLVES C G = R
00275         IF(CROUT.OR.GSEB.OR.PREBE) THEN
00276           CALL DOWNUP( G, AUX , R , 'D' , MESH )
00277           IF(NCSIZE.GT.1) CALL PARMOY(G,MESH)
00278         ELSEIF(PRE3D) THEN
00279           CALL CPSTVC(R%ADR(1)%P,G%ADR(1)%P)
00280           CALL TRID3D(AUX%X%R,G%ADR(1)%P%R,R%ADR(1)%P%R,
00281      &                MESH%NPOIN,BIEF_NBPTS(11,MESH))
00282         ENDIF
00283 !       COMPUTES RMGM AND RM1GM1
00284         RM1GM1 = STO1
00285         RMGM = P_DOTS(R,G,MESH)
00286         STO1=RMGM
00287 !
00288       ENDIF
00289 !
00290 !-----------------------------------------------------------------------
00291 ! COMPUTES D BY RECURRENCE:
00292 !-----------------------------------------------------------------------
00293 !
00294       BETA = RMGM / RM1GM1
00295 !     CALL OS( 'X=CX    ' , D , D , D , BETA )
00296 !     CALL OS( 'X=X+Y   ' , D , G , G , C    )
00297 !     OPTIMISED EQUIVALENT OF THE 2 OS CALLS:
00298       CALL OS( 'X=Y+CZ  ' , X=D , Y=G , Z=D , C=BETA )
00299 !
00300 !-----------------------------------------------------------------------
00301 ! PRECONDITIONING :
00302 !-----------------------------------------------------------------------
00303 !
00304       IF(PREC) THEN
00305 !       COMPUTES RMDM
00306         RMDM=P_DOTS(R,D,MESH)
00307       ENDIF
00308 !
00309 !-----------------------------------------------------------------------
00310 ! COMPUTES A D :
00311 !-----------------------------------------------------------------------
00312 !
00313       CALL MATRBL( 'X=AY    ',AD,A,D,C,MESH)
00314 !
00315 !-----------------------------------------------------------------------
00316 ! COMPUTES RO
00317 !-----------------------------------------------------------------------
00318 !
00319       DAD = P_DOTS(D,AD,MESH)
00320       RO = RMDM/DAD
00321 !
00322 ! COMPUTES X(M) = X(M-1) - RO * D
00323 !
00324       CALL OS('X=X+CY  ',X=X,Y=D,C=-RO)
00325 !
00326       IF(M.LT.CFG%NITMAX) GO TO 2
00327 !
00328 !-----------------------------------------------------------------------
00329 !
00330 !     IF(INFOGR) THEN
00331         TESTL = SQRT( RMRM / XL )
00332         IF (RELAT) THEN
00333            IF (LNG.EQ.1) WRITE(LU,103) M,TESTL
00334            IF (LNG.EQ.2) WRITE(LU,104) M,TESTL
00335         ELSE
00336            IF (LNG.EQ.1) WRITE(LU,203) M,TESTL
00337            IF (LNG.EQ.2) WRITE(LU,204) M,TESTL
00338         ENDIF
00339 !     ENDIF
00340       GO TO 1000
00341 !
00342 !-----------------------------------------------------------------------
00343 !
00344 900   CONTINUE
00345 !
00346       IF(INFOGR) THEN
00347         TESTL = SQRT( RMRM / XL )
00348         IF (RELAT) THEN
00349            IF (LNG.EQ.1) WRITE(LU,101) M,TESTL
00350            IF (LNG.EQ.2) WRITE(LU,102) M,TESTL
00351         ELSE
00352            IF (LNG.EQ.1) WRITE(LU,201) M,TESTL
00353            IF (LNG.EQ.2) WRITE(LU,202) M,TESTL
00354         ENDIF
00355       ENDIF
00356 !
00357 1000  CONTINUE
00358       RETURN
00359 !
00360 !-----------------------------------------------------------------------
00361 !
00362 !   FORMATS
00363 !
00364 101   FORMAT(1X,'GRACJG (BIEF) : ',
00365      &                     1I8,' ITERATIONS, PRECISION RELATIVE:',G16.7)
00366 102   FORMAT(1X,'GRACJG (BIEF) : ',
00367      &                     1I8,' ITERATIONS, RELATIVE PRECISION:',G16.7)
00368 201   FORMAT(1X,'GRACJG (BIEF) : ',
00369      &                     1I8,' ITERATIONS, PRECISION ABSOLUE :',G16.7)
00370 202   FORMAT(1X,'GRACJG (BIEF) : ',
00371      &                     1I8,' ITERATIONS, ABSOLUTE PRECISION:',G16.7)
00372 103   FORMAT(1X,'GRACJG (BIEF) : MAX D''ITERATIONS ATTEINT:',
00373      &                     1I8,' PRECISION RELATIVE:',G16.7)
00374 104   FORMAT(1X,'GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00375      &                     1I8,' RELATIVE PRECISION:',G16.7)
00376 105   FORMAT(1X,'GRACJG (BIEF) : ',
00377      &         ' SOLUTION X=0 CAR NORME L2 DE B QUASI-NULLE:',G16.7)
00378 106   FORMAT(1X,'GRACJG (BIEF) : ',
00379      &        ' SOLUTION X=0 BECAUSE L2-NORM OF B VERY SMALL:',G16.7)
00380 203   FORMAT(1X,'GRACJG (BIEF) : MAX D''ITERATIONS ATTEINT:',
00381      &                     1I8,' PRECISION ABSOLUE :',G16.7)
00382 204   FORMAT(1X,'GRACJG (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00383      &                     1I8,' ABSOLUTE PRECISION:',G16.7)
00384 !
00385 !-----------------------------------------------------------------------
00386 !
00387       END

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