rescjg.f

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

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