1 : C BCGG SOURCE GOUNAND 11/10/07 21:15:01 7149
2 : SUBROUTINE BCGG(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ RTLD,XTLD,UHAT,R,U,Z,ZM1,Y,YP,BP,TRAV,
6 : $ ITER,INMV,BRTOL,LBCG,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
7 : IMPLICIT REAL*8 (A-H,O-Z)
8 : IMPLICIT INTEGER(I-N)
9 : C***********************************************************************
10 : C NOM : BCGG
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,ZM1.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,Y.IZA,YP.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 RHOTOL,OMETOL
98 : REAL*8 GDOT,GNRM2,PETIT,ZMXRNR,ZMXRNX
99 : LOGICAL LCPRES,LUPDAX
100 : * ..
101 : * .. External subroutines and functions..
102 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
103 : * ..
104 : * .. Executable Statements ..
105 : *
106 : *
107 : IRET = 0
108 : INMV = 0
109 : IWARN = 0
110 : IGRAN = 0
111 : MGRAN = 3
112 : MAXIT = ITER
113 : TOL = RESID
114 : *
115 : * Same memory space
116 : RHAT=UHAT
117 : *
118 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgg.eso'
119 : *
120 : * Set parameter tolerances.
121 : *
122 : RHOTOL = BRTOL
123 : SIGTOL = BRTOL
124 : RRMAX = 1.D0/(XZPREC**0.25D0)
125 : * XZPREC* 1.D0 plante si tous les termes de la matrice Z sont égaux
126 : * sur IBM
127 : PETIT = XZPREC*10.D0
128 : *
129 : * Set initial residual.
130 : *
131 : C r(1)=b
132 : RI=R.WRK(1)
132 : CALL GCOPY(B,RI)
133 : C r(1)=b-Ax
134 : * IF (GNRM2(X).NE.ZERO) THEN
135 : C r(1)=b-Ax
135 : * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
136 : * ENDIF
136 : CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,RI)
137 : INMV=INMV+1
137 : CALL GCOPY(RI,BP)
138 : *
139 : ITER = 0
140 : BNRM2 = GNRM2(B)
141 : IF(IMPR.GT.6) THEN
142 : WRITE(IOIMP,*) '||B||=',BNRM2
143 : ENDIF
144 : RNRM2 = GNRM2(RI)
145 : IF (ICALRS.EQ.1) BNRM2=RNRM2
146 : *
147 : IF (BNRM2.LT.XPETIT) BNRM2 = ONE
148 : RESID = RNRM2 / BNRM2
149 : IF(IMPR.GT.6) THEN
150 : WRITE(IOIMP,*) 'Résidu initial=',RESID
151 : ENDIF
152 : RESMAX = RESID * RRMAX
153 : LRES.PROG(ITER+1)=RESID
154 : LNMV.LECT(ITER+1)=INMV
155 : IF (RESID.LE.TOL) GOTO 30
156 : C
157 : C Variables for reliable updating
158 : C
159 : ZMXRNR=RNRM2
160 : ZMXRNX=RNRM2
161 : LCPRES=.FALSE.
162 : LUPDAX=.FALSE.
163 : XDELTA=1.D-2
164 : C r[tilde]=r(1)
164 : CALL GCOPY(RI,RTLD)
165 : C alpha = rho0 = omega = 1
166 : ALPHA = ONE
167 : RHO0 = ONE
168 : OMEGA = ONE
169 : *
170 : 10 CONTINUE
171 : *
172 : * Perform BiConjugate Gradient part.
173 : *
174 : ITER = ITER + 1
175 : IF(IMPR.GT.7) THEN
176 : WRITE(IOIMP,*) 'ITER=',ITER
177 : ENDIF
178 : *
179 : C rho0 = -omega rho0
180 : RHO0 = -OMEGA * RHO0
181 : DO J = 1 , LBCG
182 : RI = R.WRK(J)
183 : RHO1 = GDOT(RI,RTLD,IDDOT)
184 : IF (ABS(RHO0).LT.RHOTOL) THEN
185 : RHO0=SIGN(XZPREC,RHO0)
186 : IWARN=IWARN+1
187 : ENDIF
188 : BETA = ALPHA * (RHO1 / RHO0)
189 : RHO0 = RHO1
190 : DO I = 1 , J
191 : UI = U.WRK(I)
192 : RI = R.WRK(I)
192 : CALL GSCAL(-BETA,UI)
192 : CALL GAXPY(ONE,RI,UI)
193 : ENDDO
194 : UI = U.WRK(J)
195 : UIP = U.WRK(J+1)
195 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
196 : $ ,UHAT,UI)
196 : CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,UHAT,ZERO,UIP)
197 : INMV=INMV+1
198 : C
199 : SIGMA = GDOT(UIP,RTLD,IDDOT)
200 : IF (ABS(SIGMA).LT.SIGTOL) THEN
201 : SIGMA=SIGN(XZPREC,SIGMA)
202 : IWARN=IWARN+1
203 : ENDIF
204 : ALPHA = RHO1 / SIGMA
205 : UI = U.WRK(1)
205 : CALL GAXPY( ALPHA, UI , XTLD)
206 : DO I = 1 , J
207 : UIP = U.WRK(I+1)
208 : RI = R.WRK(I)
208 : CALL GAXPY(-ALPHA,UIP,RI)
209 : ENDDO
210 : RI = R.WRK(J)
211 : RIP = R.WRK(J+1)
211 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
212 : $ ,RHAT,RI)
212 : CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,RHAT,ZERO,RIP)
213 : INMV=INMV+1
214 : ENDDO
215 : *
216 : * Perform polynomial part.
217 : *
218 : ZMAX=0.D0
219 : DO I=1,LBCG
220 : DO J=1,I
221 : RI = R.WRK(I+1)
222 : RIP = R.WRK(J+1)
223 : TAU = GDOT (RIP, RI,IDDOT)
224 : Z.IJ(I,J) = TAU
225 : IF (I .NE. J) THEN
226 : Z.IJ(J,I) = TAU
227 : ENDIF
228 : ENDDO
229 : ZMAX=MAX(ZMAX,ABS(Z.IJ(I,I)))
230 : RI = R.WRK(1)
231 : RIP = R.WRK(I+1)
232 : Y.A(I)=GDOT(RI,RIP,IDDOT)
233 : ENDDO
234 : * Ca arrive dans invdiag et invide.dgibi !
235 : IF (ZMAX.LT.(XPETIT**0.5D0)) THEN
236 : ZMAX=XZPREC
237 : ENDIF
238 : * Petite pénalisation des familles pour l'inversibilité
239 : DO I=1,LBCG
240 : Z.IJ(I,I)=Z.IJ(I,I)+PETIT*ZMAX
241 : ENDDO
242 : IIMPR=1
242 : CALL IVMAT(LBCG,Z.IJ,TRAV.LECT,ZM1.IJ,XDETZ,IIMPR,IRET)
243 : * IF (IRET.NE.0) GOTO 25
244 : IF (IRET.NE.0) THEN
245 : * SEGPRT,Z
246 : GOTO 25
247 : ENDIF
248 : *
249 : DO I=1,LBCG
250 : YP.A(I)=ZERO
251 : DO J=1,LBCG
252 : YP.A(I)=YP.A(I)+(ZM1.IJ(I,J)*Y.A(J))
253 : ENDDO
254 : ENDDO
255 : OMEGA=YP.A(LBCG)
256 : DO I=1,LBCG
257 : UI =U.WRK(1)
258 : UIP=U.WRK(I+1)
258 : CALL GAXPY (-YP.A(I),UIP,UI)
259 : RI=R.WRK(I)
259 : CALL GAXPY ( YP.A(I),RI,XTLD)
260 : RI =R.WRK(1)
261 : RIP=R.WRK(I+1)
261 : CALL GAXPY (-YP.A(I),RIP,RI)
262 : ENDDO
263 : *
264 : * Compute residual
265 : *
266 : RI=R.WRK(1)
267 : RNRM2I=GNRM2(RI)
268 : RESID=RNRM2I / BNRM2
269 : IF (RESID.GT.RESMAX) THEN
270 : IGRAN=IGRAN+1
271 : ELSE
272 : IGRAN=0
273 : ENDIF
274 : LRES.PROG(ITER+1)=RESID
275 : LNMV.LECT(ITER+1)=INMV
276 : * WRITE(IOIMP,*) 'ITER=',ITER,' RESID=',RESID
277 : *
278 : * Reliable update
279 : *
280 : ZMXRNX=MAX(RNRM2I,ZMXRNX)
281 : ZMXRNR=MAX(RNRM2I,ZMXRNR)
282 : LUPDAX=(RNRM2I.LT.XDELTA*RNRM2.AND.RNRM2.LT.ZMXRNX)
283 : LCPRES=((RNRM2I.LT.XDELTA*ZMXRNR.AND.RNRM2.LT.ZMXRNR)
284 : $ .OR.LUPDAX)
285 : IF (LCPRES) THEN
286 : RI=R.WRK(1)
286 : CALL GCOPY(BP,RI)
286 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
287 : $ ,UHAT,XTLD)
287 : CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,UHAT,ONE,RI)
288 : INMV=INMV+1
289 : ZMXRNR=RNRM2I
290 : *! WRITE(IOIMP,*) 'Compute residual'
291 : IF (LUPDAX) THEN
291 : CALL GAXPY(ONE,UHAT,X)
291 : CALL GSCAL(ZERO,XTLD)
291 : CALL GCOPY(RI,BP)
292 : ZMXRNX=RNRM2I
293 : *! WRITE(IOIMP,*) 'Update approx'
294 : ENDIF
295 : ENDIF
296 : *
297 : * Check for tolerance.
298 : *
299 : IF (RESID.LE.TOL) THEN
299 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
300 : $ ,UHAT,XTLD)
300 : CALL GAXPY(ONE,UHAT,X)
301 : GOTO 30
302 : ENDIF
303 : *
304 : * Iteration fails.
305 : *
306 : * IF (ITER.EQ.MAXIT) THEN
307 : IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
307 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
308 : $ ,UHAT,XTLD)
308 : CALL GAXPY(ONE,UHAT,X)
309 : IRET = 1
310 : IF (IMPR.GT.0) THEN
311 : C WRITE(IOIMP,*)
312 : C $ 'bcgg.eso : Convergence to tolerance not achieved'
313 : C WRITE(IOIMP,*) 'INMV=',INMV,
314 : C $ ' TOL=',TOL,
315 : C $ ' RESID=',RESID
316 : IF (IWARN.GT.0) THEN
317 : WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
318 : $ ,' times.'
319 : ENDIF
320 : ENDIF
321 : RETURN
322 : ENDIF
323 : *
324 : * Do some more
325 : *
326 : GOTO 10
327 : *
328 : * A breakdown in the method has occured.
329 : *
330 : 25 CONTINUE
331 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
332 : WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
333 : $ ,' times.'
334 : ENDIF
335 : IRET=-1
336 : RETURN
337 : *
338 : * Iteration successful
339 : *
340 : 30 CONTINUE
341 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
342 : WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
343 : $ ,' times.'
344 : ENDIF
345 : IRET = 0
346 : RETURN
347 : *
348 : * End of BCGG.
349 : *
350 : END
351 :
352 :
353 :
354 :
355 :
356 :
357 :
358 :
359 :
360 :
361 :
362 :
363 :
© Cast3M 2003 - All rights reserved.
Disclaimer