cgsqua.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\cgsqua.f
00002 !
00144                      SUBROUTINE CGSQUA
00145 !                    *****************
00146 !
00147      &(X,A,B,MESH, G,G0,P,K,H,AHPK,CFG,INFOGR)
00148 !
00149 !***********************************************************************
00150 ! BIEF   V6P1                                   21/08/2010
00151 !***********************************************************************
00152 !
00153 !
00154 !
00155 !
00156 !
00157 !
00158 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00159 !| A              |-->| MATRIX OF THE SYSTEM
00160 !| AHPK           |<->| WORK STRUCTURE
00161 !| B              |-->| RIGHT-HAND SIDE OF SYSTEM
00162 !| CFG            |-->| CFG(1): STORAGE OF MATRIX
00163 !|                |   | CFG(2): MATRIX VECTOR PRODUCT
00164 !| G              |<->| GRADIENT.
00165 !| G0             |<->| INITIAL GRADIENT
00166 !| H              |<->| WORK STRUCTURE
00167 !| INFOGR         |-->| IF YES, INFORMATION PRINTED
00168 !| K              |<->| WORK STRUCTURE
00169 !| MESH           |-->| MESH STRUCTURE
00170 !| P              |<->| WORK STRUCTURE
00171 !| X              |<--| INITIAL VALUE, THEN SOLUTION
00172 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00173 !
00174       USE BIEF, EX_CGSQUA => CGSQUA
00175 !
00176       IMPLICIT NONE
00177       INTEGER LNG,LU
00178       COMMON/INFO/LNG,LU
00179 !
00180 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00181 !
00182       TYPE(BIEF_OBJ)  , INTENT(INOUT) :: X,G,G0,P,K,H,AHPK
00183       TYPE(BIEF_OBJ)  , INTENT(IN   ) :: B
00184       TYPE(BIEF_OBJ)  , INTENT(IN)    :: A
00185       TYPE(SLVCFG)    , INTENT(INOUT) :: CFG
00186       LOGICAL         , INTENT(IN)    :: INFOGR
00187       TYPE(BIEF_MESH) , INTENT(INOUT) :: MESH
00188 !
00189 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00190 !
00191       DOUBLE PRECISION XL,RO,TESTL,RL,GMP1G0,BETA,GMG0,C
00192 !
00193       INTEGER M
00194 !
00195       LOGICAL RELAT
00196 !
00197       INTRINSIC SQRT
00198 !
00199 !-----------------------------------------------------------------------
00200 !
00201 ! COMPUTES THE NORM OF THE SECOND MEMBER
00202 !
00203       XL = P_DOTS(B,B,MESH)
00204 !
00205       IF(XL.LT.1.D0) THEN
00206         XL = 1.D0
00207         RELAT = .FALSE.
00208       ELSE
00209         RELAT = .TRUE.
00210       ENDIF
00211 !
00212       M = 0
00213 !
00214 ! INITIALISES G  : A X0 - B
00215 !
00216       CALL MATRBL( 'X=AY    ',G,A,X,C,  MESH)
00217 !
00218       CALL OS( 'X=X-Y   ' , G , B , B , C )
00219 !
00220 ! CHECKS THAT THE ACCURACY HAS NOT ALREADY BEEN REACHED
00221 !
00222       RL   = P_DOTS(G,G,MESH)
00223 !
00224       IF(RL.LT.CFG%EPS**2*XL) THEN
00225         TESTL = SQRT(RL/XL)
00226         IF (INFOGR) THEN
00227           IF(RELAT) THEN
00228             IF (LNG.EQ.1) WRITE(LU,100) M,TESTL
00229             IF (LNG.EQ.2) WRITE(LU,101) M,TESTL
00230           ELSE
00231             IF (LNG.EQ.1) WRITE(LU,200) M,TESTL
00232             IF (LNG.EQ.2) WRITE(LU,201) M,TESTL
00233           ENDIF
00234         ENDIF
00235         GOTO 1000
00236       ENDIF
00237 !
00238 ! INITIALISES G0 , P , AND K
00239 !
00240       CALL OS('X=Y     ' , G0 , G , G , C )
00241       CALL OS('X=Y     ' , P  , G , G , C )
00242       CALL OS('X=Y     ' , K  , G , G , C )
00243 !
00244       M = 1
00245 20    CONTINUE
00246 !
00247 ! COMPUTES AP (IN H, RECOMPUTED LATER)
00248 !
00249       CALL MATRBL( 'X=AY    ',H,A,P,C,  MESH)
00250 !
00251 ! COMPUTES RO
00252 !
00253       RO = P_DOTS(K,G0,MESH) / P_DOTS(H,G0,MESH)
00254 !
00255 ! COMPUTES H+K (IN H, AP ALREADY BEING IN H)
00256 !
00257       CALL OS( 'X=CX    ' , H , H , H , -RO   )
00258       CALL OS( 'X=X+CY  ' , H , K , K , 2.D0  )
00259 !
00260 !            M+1   M       M       M
00261 ! COMPUTES  X   = X  -   RO * (H+K)     (H+K IN H)
00262 !
00263       CALL OS( 'X=X+CY  ' , X , H , H , -RO )
00264 !
00265 ! COMPUTES A(H+K)  (HERE H+K IN H, A(H+K) IN AHPK)
00266 !
00267       CALL MATRBL( 'X=AY    ',AHPK,A,H,C,  MESH)
00268 !
00269 !             M     0
00270 ! COMPUTES ( G   , G  )
00271 !
00272       GMG0 = P_DOTS(G,G0,MESH)
00273 !
00274 ! COMPUTES GM
00275 !
00276       CALL OS( 'X=X+CY  ' , G , AHPK , AHPK , -RO )
00277 !
00278       RL   = P_DOTS(G,G,MESH)
00279       IF (RL.GT.CFG%EPS**2*XL) THEN
00280         IF (M.GE.CFG%NITMAX) THEN
00281             TESTL=SQRT(RL/XL)
00282             IF(RELAT) THEN
00283               IF (LNG.EQ.1) WRITE(LU,102) M,TESTL
00284               IF (LNG.EQ.2) WRITE(LU,103) M,TESTL
00285             ELSE
00286               IF (LNG.EQ.1) WRITE(LU,202) M,TESTL
00287               IF (LNG.EQ.2) WRITE(LU,203) M,TESTL
00288             ENDIF
00289           GOTO 1000
00290         ELSE
00291           M = M + 1
00292         ENDIF
00293       ELSE
00294         IF(INFOGR) THEN
00295           TESTL=SQRT(RL/XL)
00296           IF(RELAT) THEN
00297             IF (LNG.EQ.1) WRITE(LU,100) M,TESTL
00298             IF (LNG.EQ.2) WRITE(LU,101) M,TESTL
00299           ELSE
00300             IF (LNG.EQ.1) WRITE(LU,200) M,TESTL
00301             IF (LNG.EQ.2) WRITE(LU,201) M,TESTL
00302           ENDIF
00303         ENDIF
00304         GOTO 1000
00305       ENDIF
00306 !
00307 !             M+1   0
00308 ! COMPUTES ( G   , G  )
00309 !
00310       GMP1G0 = P_DOTS(G,G0,MESH)
00311 !
00312 ! COMPUTES BETA
00313 !
00314       BETA = GMP1G0 / GMG0
00315 !
00316 !           M
00317 ! COMPUTES H   (H+K IN H)
00318 !
00319       CALL OS('X=X-Y   ' , H , K , K , C )
00320 !
00321 !           M+1
00322 ! COMPUTES P
00323 !
00324       CALL OS('X=CX    ' , P , P , P , BETA**2 )
00325       CALL OS('X=X+Y   ' , P , G , G , C       )
00326       CALL OS('X=X+CY  ' , P , H , H , 2*BETA  )
00327 !
00328 !           M+1
00329 ! COMPUTES K
00330 !
00331       CALL OS('X=Y     ' , K , G , G , C    )
00332       CALL OS('X=X+CY  ' , K , H , H , BETA )
00333 !
00334       GOTO 20
00335 !
00336 1000  RETURN
00337 !
00338 !-----------------------------------------------------------------------
00339 !
00340 !   FORMATS
00341 !
00342 100   FORMAT(1X,'CGSQUA (BIEF) : ',1I8,' ITERATIONS',
00343      &          ' PRECISION RELATIVE: ',G16.7)
00344 101   FORMAT(1X,'CGSQUA (BIEF) : ',1I8,' ITERATIONS',
00345      &          ' RELATIVE PRECISION: ',G16.7)
00346 200   FORMAT(1X,'CGSQUA (BIEF) : ',1I8,' ITERATIONS',
00347      &          ' PRECISION ABSOLUE: ',G16.7)
00348 201   FORMAT(1X,'CGSQUA (BIEF) : ',1I8,' ITERATIONS',
00349      &          ' ABSOLUTE PRECISION: ',G16.7)
00350 102   FORMAT(1X,'CGSQUA (BIEF) : MAX D''ITERATIONS ATTEINT ',1I8,
00351      &          ' PRECISION RELATIVE: ',G16.7)
00352 103   FORMAT(1X,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1I8,
00353      &          ' RELATIVE PRECISION: ',G16.7)
00354 202   FORMAT(1X,'CGSQUA (BIEF) : MAX D''ITERATIONS ATTEINT ',1I8,
00355      &          ' PRECISION ABSOLUE: ',G16.7)
00356 203   FORMAT(1X,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1I8,
00357      &          ' ABSOLUTE PRECISION:',G16.7)
00358 !
00359 !-----------------------------------------------------------------------
00360 !
00361       END

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