equnor.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\equnor.f
00002 !
00084                      SUBROUTINE EQUNOR
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: A MULTIPLIED BY DESCENT GRADIENT
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_EQUNOR => EQUNOR
00114 !
00115       IMPLICIT NONE
00116       INTEGER LNG,LU
00117       COMMON/INFO/LNG,LU
00118 !
00119 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00120 !
00121       TYPE(SLVCFG), INTENT(INOUT)    :: CFG
00122       TYPE(BIEF_OBJ), INTENT(INOUT)  :: B
00123       TYPE(BIEF_OBJ), INTENT(INOUT)  :: D,AD,G,AG,R,X
00124       TYPE(BIEF_MESH), INTENT(INOUT) :: MESH
00125       TYPE(BIEF_OBJ), INTENT(IN)     :: A
00126       TYPE(BIEF_OBJ), INTENT(INOUT)  :: AUX
00127       LOGICAL, INTENT(IN)            :: INFOGR
00128 !
00129 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00130 !
00131       INTEGER M
00132 !
00133       DOUBLE PRECISION XL,RMRM,TESTL
00134       DOUBLE PRECISION BETA,ADAD,RO
00135       DOUBLE PRECISION TGMTGM,TG1TG1,C
00136 !
00137       LOGICAL RELAT,PREC,CROUT,GSEB
00138 !
00139 !-----------------------------------------------------------------------
00140 !
00141       INTRINSIC SQRT
00142 !
00143 !-----------------------------------------------------------------------
00144 !
00145 !   INITIALISES
00146 !
00147       CROUT =.FALSE.
00148       IF(7*(CFG%PRECON/7).EQ.CFG%PRECON) CROUT=.TRUE.
00149       GSEB=.FALSE.
00150       IF(11*(CFG%PRECON/11).EQ.CFG%PRECON) GSEB=.TRUE.
00151       PREC=.FALSE.
00152       IF(CROUT.OR.GSEB.OR.13*(CFG%PRECON/13).EQ.CFG%PRECON) PREC=.TRUE.
00153 !
00154 !-----------------------------------------------------------------------
00155 !   INITIALISES
00156 !-----------------------------------------------------------------------
00157 !
00158       M   = 0
00159 !
00160 !  NORMALISES THE SECOND MEMBER TO COMPUTE THE RELATIVE PRECISION:
00161 !
00162       XL = P_DOTS(B,B,MESH)
00163 !
00164       IF(XL.LT.1.D0) THEN
00165         XL = 1.D0
00166         RELAT = .FALSE.
00167       ELSE
00168         RELAT = .TRUE.
00169       ENDIF
00170 !
00171 ! COMPUTES THE INITIAL RESIDUAL AND POSSIBLY EXITS:
00172 !
00173       CALL MATRBL( 'X=AY    ',R,A,X,  C,MESH)
00174 !
00175       CALL OS( 'X=X-Y   ' , R , B , B , C )
00176       RMRM   = P_DOTS(R,R,MESH)
00177 !
00178       IF (RMRM.LT.CFG%EPS**2*XL) GO TO 900
00179 !
00180 !-----------------------------------------------------------------------
00181 ! PRECONDITIONING :
00182 !-----------------------------------------------------------------------
00183 !
00184       IF(PREC) THEN
00185 !
00186 !       COMPUTES C G0 = R
00187         CALL DOWNUP(G, AUX , R , 'D' , MESH)
00188 !
00189 !  C IS HERE CONSIDERED SYMMETRICAL,
00190 !  SHOULD OTHERWISE SOLVE TC GPRIM = G
00191 !
00192 !        T -1
00193 !         C   G   IS PUT IN B
00194 !
00195         CALL DOWNUP(B , AUX , G , 'T' , MESH)
00196 !
00197       ENDIF
00198 !
00199 !-----------------------------------------------------------------------
00200 ! COMPUTES THE DIRECTION OF INITIAL DESCENT:
00201 !-----------------------------------------------------------------------
00202 !
00203       IF(PREC) THEN
00204         CALL MATRBL( 'X=TAY   ',D,A,B,  C,MESH)
00205       ELSE
00206         CALL MATRBL( 'X=TAY   ',D,A,G,  C,MESH)
00207       ENDIF
00208 !
00209       TGMTGM = P_DOTS(D,D,MESH)
00210 !
00211 !-----------------------------------------------------------------------
00212 ! COMPUTES THE INITIAL PRODUCT A D:
00213 !-----------------------------------------------------------------------
00214 !
00215       CALL MATRBL('X=AY    ',AD,A,D,C,MESH)
00216 !
00217 !-----------------------------------------------------------------------
00218 !
00219       IF(PREC) THEN
00220 !
00221 !         COMPUTES  C DPRIM = AD  (DPRIM PUT IN AG)
00222           CALL DOWNUP(AG, AUX , AD , 'D' , MESH)
00223 !
00224       ENDIF
00225 !
00226 !-----------------------------------------------------------------------
00227 ! COMPUTES INITIAL RO :
00228 !-----------------------------------------------------------------------
00229 !
00230       IF(PREC) THEN
00231         ADAD = P_DOTS(AG,AG,MESH)
00232       ELSE
00233         ADAD = P_DOTS(AD,AD,MESH)
00234       ENDIF
00235       RO = TGMTGM/ADAD
00236 !
00237 !-----------------------------------------------------------------------
00238 !
00239 ! COMPUTES X1 = X0 - RO  * D
00240 !
00241       CALL OS( 'X=X+CY  ' , X , D , D , -RO )
00242 !
00243 !-----------------------------------------------------------------------
00244 !  ITERATIONS LOOP:
00245 !-----------------------------------------------------------------------
00246 !
00247 2     M  = M  + 1
00248 !
00249 !-----------------------------------------------------------------------
00250 ! COMPUTES THE RESIDUAL : R(M) = R(M-1) - RO(M-1) A D(M-1)
00251 !-----------------------------------------------------------------------
00252 !
00253       CALL OS( 'X=X+CY  ' , R , AD , AD , -RO )
00254 !
00255 !  SOME VALUES WILL CHANGE IN CASE OF PRECONDITIONING
00256 !
00257       RMRM   = P_DOTS(R,R,MESH)
00258 !
00259 ! CHECKS END :
00260 !
00261       IF (RMRM.LE.XL*CFG%EPS**2) GO TO 900
00262 !
00263 !-----------------------------------------------------------------------
00264 ! PRECONDITIONING : SOLVES C G = R
00265 !-----------------------------------------------------------------------
00266 !
00267       IF(PREC) THEN
00268 !
00269 !         UPDATES G BY RECURRENCE (IN AG: DPRIM)
00270           CALL OS( 'X=X+CY  ' , G , AG , AG , -RO )
00271 !
00272           CALL DOWNUP(B , AUX , G , 'T' , MESH)
00273 !
00274       ENDIF
00275 !
00276 !-----------------------------------------------------------------------
00277 ! COMPUTES D BY RECURRENCE:
00278 !-----------------------------------------------------------------------
00279 !
00280       IF(PREC) THEN
00281 !                                          T  T -1          T -1
00282 !                               AD IS HERE  A  C  G    B IS  C   G
00283         CALL MATRBL( 'X=TAY   ',AD,A,B,  C,MESH)
00284       ELSE
00285 !                               AD IS HERE TAG
00286         CALL MATRBL( 'X=TAY   ',AD,A,G,  C,MESH)
00287       ENDIF
00288 !
00289       TG1TG1 = TGMTGM
00290       TGMTGM = P_DOTS(AD,AD,MESH)
00291       BETA = TGMTGM / TG1TG1
00292 !
00293       CALL OS( 'X=CX    ' , D , D , D , BETA )
00294 !
00295 !                               AD IS HERE TAG
00296       CALL OS( 'X=X+Y   ' , D , AD , AD , C   )
00297 !
00298 !-----------------------------------------------------------------------
00299 ! COMPUTES A D :
00300 !-----------------------------------------------------------------------
00301 !
00302       CALL MATRBL( 'X=AY    ',AD,A,D,  C,MESH)
00303 !
00304       IF(PREC) THEN
00305 !
00306 !           COMPUTES  C DPRIM = AD  (DPRIM PUT IN AG)
00307             CALL DOWNUP(AG , AUX , AD , 'D' , MESH)
00308 !
00309       ENDIF
00310 !
00311 !-----------------------------------------------------------------------
00312 ! COMPUTES RO
00313 !-----------------------------------------------------------------------
00314 !
00315       IF(PREC) THEN
00316         ADAD = P_DOTS(AG,AG,MESH)
00317       ELSE
00318         ADAD = P_DOTS(AD,AD,MESH)
00319       ENDIF
00320       RO = TGMTGM/ADAD
00321 !
00322 ! COMPUTES X(M) = X(M-1) - RO * D
00323 !
00324       CALL OS( 'X=X+CY  ' , X , D , D , -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  RETURN
00358 !
00359 !-----------------------------------------------------------------------
00360 !
00361 !   FORMATS
00362 !
00363 101   FORMAT(1X,'EQUNOR (BIEF) : ',
00364      &                     1I8,' ITERATIONS, PRECISION RELATIVE:',G16.7)
00365 102   FORMAT(1X,'EQUNOR (BIEF) : ',
00366      &                     1I8,' ITERATIONS, RELATIVE PRECISION:',G16.7)
00367 201   FORMAT(1X,'EQUNOR (BIEF) : ',
00368      &                     1I8,' ITERATIONS, PRECISION ABSOLUE :',G16.7)
00369 202   FORMAT(1X,'EQUNOR (BIEF) : ',
00370      &                     1I8,' ITERATIONS, ABSOLUTE PRECISION:',G16.7)
00371 103   FORMAT(1X,'EQUNOR (BIEF) : MAX D'' ITERATIONS ATTEINT:',
00372      &                     1I8,' PRECISION RELATIVE:',G16.7)
00373 104   FORMAT(1X,'EQUNOR (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00374      &                     1I8,' RELATIVE PRECISION:',G16.7)
00375 203   FORMAT(1X,'EQUNOR (BIEF) : MAX D'' ITERATIONS ATTEINT:',
00376      &                     1I8,' PRECISION ABSOLUE :',G16.7)
00377 204   FORMAT(1X,'EQUNOR (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
00378      &                     1I8,' ABSOLUTE PRECISON:',G16.7)
00379 !
00380 !-----------------------------------------------------------------------
00381 !
00382       END

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