cgstab.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\cgstab.f
00002 !
00079                      SUBROUTINE CGSTAB
00080 !                    *****************
00081 !
00082      &(X, A,B , MESH, P,Q,R,S,T,V, CFG,INFOGR,AUX)
00083 !
00084 !***********************************************************************
00085 ! BIEF   V6P1                                   21/08/2010
00086 !***********************************************************************
00087 !
00088 !
00089 !
00090 !
00091 !
00092 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00093 !| A              |-->| MATRIX OF THE SYSTEM
00094 !| AUX            |-->| PRECONDITIONING MATRIX
00095 !| B              |-->| RIGHT-HAND SIDE OF SYSTEM
00096 !| CFG            |-->| CFG(1): STORAGE OF MATRIX
00097 !|                |   | CFG(2): MATRIX VECTOR PRODUCT
00098 !| INFOGR         |-->| IF YES, INFORMATION PRINTED
00099 !| MESH           |-->| MESH STRUCTURE
00100 !| P              |<->| WORK STRUCTURE
00101 !| Q              |<->| WORK STRUCTURE
00102 !| R              |<->| WORK STRUCTURE
00103 !| S              |<->| WORK STRUCTURE
00104 !| T              |<->| WORK STRUCTURE
00105 !| V              |<->| WORK STRUCTURE
00106 !| X              |<--| VALEUR INITIALE, PUIS SOLUTION
00107 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00108 !
00109       USE BIEF, EX_CGSTAB => CGSTAB
00110 !
00111       IMPLICIT NONE
00112       INTEGER LNG,LU
00113       COMMON/INFO/LNG,LU
00114 !
00115 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00116 !
00117       TYPE(BIEF_OBJ)  , INTENT(INOUT) :: X,P,Q,R,S,T,V
00118       TYPE(BIEF_OBJ)  , INTENT(IN)    :: AUX,A,B
00119       TYPE(BIEF_MESH) , INTENT(INOUT) :: MESH
00120       TYPE(SLVCFG)    , INTENT(INOUT) :: CFG
00121       LOGICAL         , INTENT(IN)    :: INFOGR
00122 !
00123 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00124 !
00125       DOUBLE PRECISION ALFA,ALFA1,BETA,BETA1,OMEG,OMEG1,OMEG2
00126       DOUBLE PRECISION XL,TESTL,RMRM,C
00127 !
00128       INTEGER M
00129 !
00130       LOGICAL RELAT,CROUT
00131 !
00132       INTRINSIC SQRT
00133 !
00134 !-----------------------------------------------------------------------
00135 !
00136 !   INITIALISES
00137       CROUT=.FALSE.
00138       IF(7*(CFG%PRECON/7).EQ.CFG%PRECON) CROUT=.TRUE.
00139 !
00140 !-----------------------------------------------------------------------
00141 !   INITIALISES
00142 !-----------------------------------------------------------------------
00143 !
00144       M   = 0
00145 !
00146 !  NORM OF THE SECOND MEMBER TO COMPUTE THE RELATIVE ACCURACY:
00147 !
00148       XL = P_DOTS(B,B,MESH)
00149       IF (XL.LT.1.D0) THEN
00150         XL = 1.D0
00151         RELAT = .FALSE.
00152       ELSE
00153         RELAT = .TRUE.
00154       ENDIF
00155 !
00156 ! IF THE SECOND MEMBER IS 0, X=0 AND IT'S THE END
00157 !
00158       IF(XL.LT.CFG%ZERO**2) THEN
00159         RMRM = 0.D0
00160         CALL OS( 'X=C     ' , X , X , X , 0.D0 )
00161         GOTO 900
00162       ENDIF
00163 !
00164 ! COMPUTES THE INITIAL RESIDUAL AND EXITS IF RESIDUAL IS SMALL:
00165 !
00166       CALL MATRBL( 'X=AY    ',V,A,X,C,  MESH)
00167 !
00168       CALL OS( 'X=Y-Z   ' , R , B , V , C )
00169       RMRM   = P_DOTS(R,R,MESH)
00170 !
00171       IF (RMRM.LT.CFG%EPS**2*XL) GO TO 900
00172 !
00173 !-----------------------------------------------------------------------
00174 ! PRECONDITIONING :
00175 !-----------------------------------------------------------------------
00176 !
00177       IF(CROUT) THEN
00178 !       COMPUTES C R  = B
00179         CALL DOWNUP(R, AUX , B , 'D' , MESH)
00180         IF(NCSIZE.GT.1) CALL PARMOY(R,MESH)
00181       ELSE
00182         CALL OS( 'X=Y     ' , R , B , B , C )
00183       ENDIF
00184 !
00185 !-----------------------------------------------------------------------
00186 ! RESUMES INITIALISATIONS
00187 !-----------------------------------------------------------------------
00188 !
00189       IF(CROUT) THEN
00190         CALL DOWNUP(V, AUX , V , 'D' , MESH)
00191         IF(NCSIZE.GT.1) CALL PARMOY(V,MESH)
00192       ENDIF
00193 !
00194       CALL OS( 'X=X-Y   ' , R , V , V , C    )
00195       CALL OS( 'X=Y     ' , P , R , R , C    )
00196       CALL OS( 'X=C     ' , V , V , V , 0.D0 )
00197       CALL OS( 'X=C     ' , Q , Q , Q , 0.D0 )
00198 !
00199       ALFA  = 1.D0
00200       BETA  = 1.D0
00201       OMEG1 = 1.D0
00202 !
00203 !-----------------------------------------------------------------------
00204 !  ITERATIONS:
00205 !-----------------------------------------------------------------------
00206 !
00207 2     M  = M  + 1
00208 !
00209       BETA1 = P_DOTS(R,P,MESH)
00210       OMEG2 = OMEG1*BETA1/BETA
00211       OMEG  = OMEG2/ALFA
00212       BETA  = BETA1
00213 !
00214       CALL OS( 'X=Y+CZ  ' , Q , R    , Q ,  OMEG )
00215       CALL OS( 'X=X+CY  ' , Q , V    , V , -OMEG2)
00216 !
00217       CALL MATRBL( 'X=AY    ',V,A,Q,C,  MESH)
00218 !
00219       IF(CROUT) THEN
00220         CALL DOWNUP(V, AUX , V , 'D' , MESH)
00221         IF(NCSIZE.GT.1) CALL PARMOY(V,MESH)
00222       ENDIF
00223 !
00224       OMEG1 = P_DOTS(P,V,MESH)
00225       OMEG1 = BETA1/OMEG1
00226 !
00227       CALL OS( 'X=Y+CZ  ' , S , R    , V , -OMEG1)
00228 !
00229       CALL MATRBL( 'X=AY    ',T,A,S,C,  MESH)
00230 !
00231       IF(CROUT) THEN
00232         CALL DOWNUP(T, AUX , T , 'D' , MESH)
00233         IF(NCSIZE.GT.1) CALL PARMOY(T,MESH)
00234       ENDIF
00235 !
00236       ALFA  = P_DOTS(T,S,MESH)
00237       ALFA1 = P_DOTS(T,T,MESH)
00238       ALFA  = ALFA/ALFA1
00239 !
00240       CALL OS( 'X=X+CY  ' , X , Q , Q ,  OMEG1)
00241       CALL OS( 'X=X+CY  ' , X , S , S ,  ALFA )
00242 !
00243       CALL OS( 'X=Y+CZ  ' , R , S , T , -ALFA )
00244 !
00245       RMRM   = P_DOTS(R,R,MESH)
00246 !
00247 ! TESTS CONVERGENCE :
00248 !
00249       IF (RMRM.LE.XL*CFG%EPS**2) GO TO 900
00250 !
00251       IF(M.LT.CFG%NITMAX) GO TO 2
00252 !
00253 !-----------------------------------------------------------------------
00254 !
00255 !     IF(INFOGR) THEN
00256         TESTL = SQRT( RMRM / XL )
00257         IF(RELAT) THEN
00258           IF (LNG.EQ.1) WRITE(LU,102) M,TESTL
00259           IF (LNG.EQ.2) WRITE(LU,104) M,TESTL
00260         ELSE
00261           IF (LNG.EQ.1) WRITE(LU,202) M,TESTL
00262           IF (LNG.EQ.2) WRITE(LU,204) M,TESTL
00263         ENDIF
00264 !     ENDIF
00265       GO TO 1000
00266 !
00267 !-----------------------------------------------------------------------
00268 !
00269 900   CONTINUE
00270 !
00271       IF(INFOGR) THEN
00272         TESTL = SQRT( RMRM / XL )
00273         IF(RELAT) THEN
00274           IF (LNG.EQ.1) WRITE(LU,101) M,TESTL
00275           IF (LNG.EQ.2) WRITE(LU,103) M,TESTL
00276         ELSE
00277           IF (LNG.EQ.1) WRITE(LU,201) M,TESTL
00278           IF (LNG.EQ.2) WRITE(LU,203) M,TESTL
00279         ENDIF
00280       ENDIF
00281 !
00282 1000  RETURN
00283 !
00284 !-----------------------------------------------------------------------
00285 !
00286 !   FORMATS
00287 !
00288 101   FORMAT(1X,'CGSTAB (BIEF) : ',1I8,' ITERATIONS',
00289      &          ' PRECISION RELATIVE: ',G16.7)
00290 103   FORMAT(1X,'CGSTAB (BIEF) : ',1I8,' ITERATIONS',
00291      &          ' RELATIVE PRECISION: ',G16.7)
00292 201   FORMAT(1X,'CGSTAB (BIEF) : ',1I8,' ITERATIONS',
00293      &          ' PRECISION ABSOLUE: ',G16.7)
00294 203   FORMAT(1X,'CGSTAB (BIEF) : ',1I8,' ITERATIONS',
00295      &          ' ABSOLUTE PRECISION: ',G16.7)
00296 102   FORMAT(1X,'CGSTAB (BIEF) : MAX D''ITERATIONS ATTEINT ',1I8,
00297      &          ' PRECISION RELATIVE: ',G16.7)
00298 104   FORMAT(1X,'CGSTAB (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1I8,
00299      &          ' RELATIVE PRECISION: ',G16.7)
00300 202   FORMAT(1X,'CGSTAB (BIEF) : MAX D''ITERATIONS ATTEINT ',1I8,
00301      &          ' PRECISION ABSOLUE: ',G16.7)
00302 204   FORMAT(1X,'CGSTAB (BIEF) : EXCEEDING MAXIMUM ITERATIONS',1I8,
00303      &          ' ABSOLUTE PRECISION:',G16.7)
00304 !
00305 !-----------------------------------------------------------------------
00306 !
00307       END

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