1 : C BCGL SOURCE GOUNAND 11/07/21 21:15:08 7046
2 : SUBROUTINE BCGL(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ RTLD,XTLD,UHAT,R,U,Z,ZM1,Y,YP,TRAV,
6 : $ ITER,INMV,BRTOL,LBCG,RESID,ICALRS,IMPR,IRET)
7 : IMPLICIT REAL*8 (A-H,O-Z)
8 : IMPLICIT INTEGER(I-N)
9 : C***********************************************************************
10 : C NOM : BCGL
11 : C DESCRIPTION :
12 : C Résolution d'un système linéaire Ax=b
13 : C par une méthode BiCGSTAB(l) préconditionnée ou non.
14 : C
15 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
16 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
17 : C mél : gounand@semt2.smts.cea.fr
18 : C REFERENCE (bibtex-like) :
19 : C @BOOK{templates,
20 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
21 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
22 : C H. Van der Vorst},
23 : C TITLE={Templates for the Solution of Linear Systems :
24 : C Building Blocks for Iterative Methods},
25 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
26 : C -> URL : http://www.netlib.org/templates/Templates.html
27 : C@TechReport{fokkema,
28 : C author = {DR Fokkema},
29 : C title = {Enhanced implementation of BiCGSTAB(l) for solving
30 : C linear systems of equations},
31 : C institution = {?},
32 : C year = {1996}}
33 : C***********************************************************************
34 : C ENTREES :
35 : C ENTREES/SORTIES :
36 : C SORTIES :
37 : C CODE RETOUR (IRET) : 0 si ok
38 : C <0 si problème (non utilisé !)
39 : C >0 si l'algorithme n'a pas convergé
40 : C***********************************************************************
41 : C VERSION : v1, 22/02/06, version initiale
42 : C HISTORIQUE : v1, 22/02/06, création
43 : C HISTORIQUE :
44 : C HISTORIQUE :
45 : C***********************************************************************
46 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
47 : C en cas de modification de ce sous-programme afin de faciliter
48 : C la maintenance !
49 : C***********************************************************************
50 : *
51 : * .. Includes .. a commenter si mise au point ok
52 : -INC CCREEL
53 : -INC CCOPTIO
54 : -INC SMLREEL
55 : POINTEUR LRES.MLREEL
56 : -INC SMLENTI
57 : POINTEUR LNMV.MLENTI
58 : POINTEUR TRAV.MLENTI
59 : -INC SMMATRIK
60 : POINTEUR KMORS.PMORS
61 : POINTEUR KISA.IZA
62 : POINTEUR INVDIA.IZA
63 : POINTEUR ILUM.PMORS
64 : POINTEUR ILUI.IZA
65 : POINTEUR B.IZA
66 : POINTEUR X.IZA
67 : * .. Work Vectors for BiCGSTAB(l)
68 : SEGMENT SPACE
69 : REAL*8 IJ(NI,NJ)
70 : ENDSEGMENT
71 : POINTEUR Z.SPACE,ZM1.SPACE
72 : SEGMENT SPACE2
73 : POINTEUR WRK(NIZA).IZA
74 : ENDSEGMENT
75 : POINTEUR R.SPACE2,U.SPACE2
76 : * .. Work Vectors for BiCGStab
77 : POINTEUR RTLD.IZA,XTLD.IZA
78 : POINTEUR RI.IZA,UI.IZA,UIP.IZA,RIP.IZA
79 : POINTEUR UHAT.IZA,RHAT.IZA,Y.IZA,YP.IZA
80 : *
81 : INTEGER IMPR
82 : * .. Scalar Arguments ..
83 : INTEGER N,NNZ
84 : INTEGER ITER,IRET
85 : REAL*8 BRTOL,RESID
86 : * ..
87 : *
88 : *
89 : * .. Variables locales
90 : * .. Parameters
91 : REAL*8 ZERO,ONE
92 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
93 : * ..
94 : INTEGER MAXIT,III
95 : REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
96 : REAL*8 RHOTOL,OMETOL
97 : REAL*8 GDOT,GNRM2,PETIT
98 : * ..
99 : * .. External subroutines and functions..
100 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
101 : * ..
102 : * .. Executable Statements ..
103 : *
104 : *
105 : IRET = 0
106 : INMV = 0
107 : IWARN = 0
108 : IGRAN = 0
109 : MGRAN = 3
110 : MAXIT = ITER
111 : TOL = RESID
112 : *
113 : * Same memory space
114 : RHAT=UHAT
115 : *
116 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgl.eso'
117 : *
118 : * Set parameter tolerances.
119 : *
120 : RHOTOL = BRTOL
121 : SIGTOL = BRTOL
122 : RRMAX = 1.D0/(XZPREC**0.25D0)
123 : PETIT = XZPREC*10.D0
124 : *
125 : * Set initial residual.
126 : *
127 : C r(1)=b
128 : RI=R.WRK(1)
128 : CALL GCOPY(B,RI)
129 : C r(1)=b-Ax
130 : * IF (GNRM2(X).NE.ZERO) THEN
131 : C r(1)=b-Ax
131 : * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
132 : * ENDIF
132 : CALL GMOMV('N',-ONE,KMORS,KISA,X,ONE,RI)
133 : INMV=INMV+1
134 : *
135 : ITER = 0
136 : BNRM2 = GNRM2(B)
137 : IF(IMPR.GT.6) THEN
138 : WRITE(IOIMP,*) '||B||=',BNRM2
139 : ENDIF
140 : RNRM2 = GNRM2(RI)
141 : IF (ICALRS.EQ.1) BNRM2=RNRM2
142 : *
143 : IF (BNRM2.LT.XPETIT) BNRM2 = ONE
144 : RESID = RNRM2 / BNRM2
145 : IF(IMPR.GT.6) THEN
146 : WRITE(IOIMP,*) 'Résidu initial=',RESID
147 : ENDIF
148 : RESMAX = RESID * RRMAX
149 : LRES.PROG(ITER+1)=RESID
150 : LNMV.LECT(ITER+1)=INMV
151 : IF (RESID.LE.TOL) GOTO 30
152 : C r[tilde]=r(1)
152 : CALL GCOPY(RI,RTLD)
153 : C alpha = rho0 = omega = 1
154 : ALPHA = ONE
155 : RHO0 = ONE
156 : OMEGA = ONE
157 : *
158 : 10 CONTINUE
159 : *
160 : * Perform BiConjugate Gradient part.
161 : *
162 : ITER = ITER + 1
163 : IF(IMPR.GT.7) THEN
164 : WRITE(IOIMP,*) 'ITER=',ITER
165 : ENDIF
166 : *
167 : C rho0 = -omega rho0
168 : RHO0 = -OMEGA * RHO0
169 : DO J = 1 , LBCG
170 : RI = R.WRK(J)
171 : RHO1 = GDOT(RI,RTLD)
172 : IF (ABS(RHO0).LT.RHOTOL) THEN
173 : RHO0=SIGN(XZPREC,RHO0)
174 : IWARN=IWARN+1
175 : ENDIF
176 : BETA = ALPHA * (RHO1 / RHO0)
177 : RHO0 = RHO1
178 : DO I = 1 , J
179 : UI = U.WRK(I)
180 : RI = R.WRK(I)
180 : CALL GSCAL(-BETA,UI)
180 : CALL GAXPY(ONE,RI,UI)
181 : ENDDO
182 : UI = U.WRK(J)
183 : UIP = U.WRK(J+1)
183 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
184 : $ ,UHAT,UI)
184 : CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,UIP)
185 : INMV=INMV+1
186 : C
187 : SIGMA = GDOT(UIP,RTLD)
188 : IF (ABS(SIGMA).LT.SIGTOL) THEN
189 : SIGMA=SIGN(XZPREC,SIGMA)
190 : IWARN=IWARN+1
191 : ENDIF
192 : ALPHA = RHO1 / SIGMA
193 : UI = U.WRK(1)
193 : CALL GAXPY( ALPHA, UI , XTLD)
194 : DO I = 1 , J
195 : UIP = U.WRK(I+1)
196 : RI = R.WRK(I)
196 : CALL GAXPY(-ALPHA,UIP,RI)
197 : ENDDO
198 : RI = R.WRK(J)
199 : RIP = R.WRK(J+1)
199 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
200 : $ ,RHAT,RI)
200 : CALL GMOMV('N',ONE,KMORS,KISA,RHAT,ZERO,RIP)
201 : INMV=INMV+1
202 : ENDDO
203 : *
204 : * Perform polynomial part.
205 : *
206 : ZMAX=0.D0
207 : DO I=1,LBCG
208 : DO J=1,I
209 : RI = R.WRK(I+1)
210 : RIP = R.WRK(J+1)
211 : TAU = GDOT (RIP, RI)
212 : Z.IJ(I,J) = TAU
213 : IF (I .NE. J) THEN
214 : Z.IJ(J,I) = TAU
215 : ENDIF
216 : ENDDO
217 : ZMAX=MAX(ZMAX,ABS(Z.IJ(I,I)))
218 : RI = R.WRK(1)
219 : RIP = R.WRK(I+1)
220 : Y.A(I)=GDOT(RI,RIP)
221 : ENDDO
222 : * Petite pénalisation des familles pour l'inversibilité
223 : DO I=1,LBCG
224 : Z.IJ(I,I)=Z.IJ(I,I)+PETIT*ZMAX
225 : ENDDO
226 : IIMPR=1
226 : CALL IVMAT(LBCG,Z.IJ,TRAV.LECT,ZM1.IJ,XDETZ,IIMPR,IRET)
227 : IF (IRET.NE.0) GOTO 25
228 : *
229 : DO I=1,LBCG
230 : YP.A(I)=ZERO
231 : DO J=1,LBCG
232 : YP.A(I)=YP.A(I)+(ZM1.IJ(I,J)*Y.A(J))
233 : ENDDO
234 : ENDDO
235 : OMEGA=YP.A(LBCG)
236 : DO I=1,LBCG
237 : UI =U.WRK(1)
238 : UIP=U.WRK(I+1)
238 : CALL GAXPY (-YP.A(I),UIP,UI)
239 : RI=R.WRK(I)
239 : CALL GAXPY ( YP.A(I),RI,XTLD)
240 : RI =R.WRK(1)
241 : RIP=R.WRK(I+1)
241 : CALL GAXPY (-YP.A(I),RIP,RI)
242 : ENDDO
243 : *
244 : * Check for tolerance.
245 : *
246 : RI=R.WRK(1)
247 : RESID=GNRM2(RI) / BNRM2
248 : IF (RESID.GT.RESMAX) THEN
249 : IGRAN=IGRAN+1
250 : ELSE
251 : IGRAN=0
252 : ENDIF
253 : LRES.PROG(ITER+1)=RESID
254 : LNMV.LECT(ITER+1)=INMV
255 : * WRITE(IOIMP,*) 'ITER=',ITER,' RESID=',RESID
256 : IF (RESID.LE.TOL) THEN
256 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
257 : $ ,UHAT,XTLD)
257 : CALL GAXPY(ONE,UHAT,X)
258 : GOTO 30
259 : ENDIF
260 : *
261 : * Iteration fails.
262 : *
263 : * IF (ITER.EQ.MAXIT) THEN
264 : IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
265 : IRET = 1
266 : IF (IMPR.GT.0) THEN
267 : C WRITE(IOIMP,*)
268 : C $ 'bcgl.eso : Convergence to tolerance not achieved'
269 : C WRITE(IOIMP,*) 'INMV=',INMV,
270 : C $ ' TOL=',TOL,
271 : C $ ' RESID=',RESID
272 : IF (IWARN.GT.0) THEN
273 : WRITE(IOIMP,*) 'BiCGStab(l) nearly failed ',IWARN
274 : $ ,' times.'
275 : ENDIF
276 : ENDIF
277 : RETURN
278 : ENDIF
279 : *
280 : * Do some more
281 : *
282 : GOTO 10
283 : *
284 : * A breakdown in the method has occured.
285 : *
286 : 25 CONTINUE
287 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
288 : WRITE(IOIMP,*) 'BiCGStab(l) nearly failed ',IWARN
289 : $ ,' times.'
290 : ENDIF
291 : IRET=-1
292 : RETURN
293 : *
294 : * Iteration successful
295 : *
296 : 30 CONTINUE
297 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
298 : WRITE(IOIMP,*) 'BiCGStab(l) nearly failed ',IWARN
299 : $ ,' times.'
300 : ENDIF
301 : IRET = 0
302 : RETURN
303 : *
304 : * End of BCGL.
305 : *
306 : END
307 :
308 :
309 :
310 :
311 :
312 :
313 :
314 :
315 :
316 :
317 :
© Cast3M 2003 - All rights reserved.
Disclaimer