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