1 : C BCG2 SOURCE GOUNAND 11/07/21 21:15:01 7046
2 : SUBROUTINE BCG2(KMORS,KISA,B,X,
3 : $ KPREC,INVDIA,ILUM,ILUI,
4 : $ LRES,LNMV,
5 : $ R,RTLD,P,PHAT,
6 : $ V,S,SHAT,T,
7 : $ ITER,INMV,BRTOL,RESID,ICALRS,IMPR,IRET)
8 : IMPLICIT REAL*8 (A-H,O-Z)
9 : IMPLICIT INTEGER(I-N)
10 : C***********************************************************************
11 : C NOM : BCG2
12 : C DESCRIPTION :
13 : C Résolution d'un système linéaire Ax=b
14 : C par la méthode du bi-gradient conjugué stabilisé
15 : C (Bi-CGSTAB) préconditionné.
16 : C
17 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
18 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
19 : C mél : gounand@semt2.smts.cea.fr
20 : C REFERENCE (bibtex-like) :
21 : C @BOOK{templates,
22 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
23 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
24 : C H. Van der Vorst},
25 : C TITLE={Templates for the Solution of Linear Systems :
26 : C Building Blocks for Iterative Methods},
27 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
28 : C -> URL : http://www.netlib.org/templates/Templates.html
29 : C***********************************************************************
30 : C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
31 : C GNRM2, GDOT (functions)
32 : C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
33 : C PSOLVE : Preconditionner solve
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, 09/02/06
43 : C HISTORIQUE : v1, 09/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 : -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
68 : POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
69 : POINTEUR V.IZA,S.IZA,SHAT.IZA,T.IZA
70 : *
71 : INTEGER IMPR
72 : * .. Scalar Arguments ..
73 : INTEGER N,NNZ
74 : INTEGER ITER,IRET
75 : REAL*8 BRTOL,RESID
76 : * ..
77 : *
78 : *
79 : * .. Variables locales
80 : * .. Parameters
81 : REAL*8 ZERO,ONE
82 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
83 : * ..
84 : INTEGER MAXIT,III
85 : REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA,TAU
86 : REAL*8 RHOTOL,OMETOL,TAUTOL
87 : REAL*8 GDOT,GNRM2
88 : * ..
89 : * .. External subroutines and functions..
90 : EXTERNAL GAXPY,GCOPY,GDOT,GNRM2,GSCAL
91 : * ..
92 : * .. Executable Statements ..
93 : *
94 : *
95 : IRET = 0
96 : INMV = 0
97 : IWARN = 0
98 : IGRAN = 0
99 : MGRAN = 5
100 : MAXIT = ITER
101 : TOL = RESID
102 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcg2.eso'
103 : *
104 : * Set parameter tolerances.
105 : *
106 : RHOTOL = BRTOL
107 : OMETOL = BRTOL
108 : TAUTOL = BRTOL
109 : RRMAX = 1.D0/(XZPREC**0.25D0)
110 : * RRMAX = 10.D0
111 : *
112 : * Set initial residual.
113 : *
114 : C r(0)=b
114 : CALL GCOPY(B,R)
115 : C r(0)=b-Ax
116 : * IF (GNRM2(X).NE.ZERO) THEN
117 : C r(0)=b-Ax
117 : * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
118 : * ENDIF
118 : CALL GMOMV('N',-ONE,KMORS,KISA,X,ONE,R)
119 : INMV=INMV+1
120 : *
121 : ITER = 0
122 : BNRM2 = GNRM2(B)
123 : IF(IMPR.GT.6) THEN
124 : WRITE(IOIMP,*) '||B||=',BNRM2
125 : ENDIF
126 : RNRM2 = GNRM2(R)
127 : IF (ICALRS.EQ.1) BNRM2=RNRM2
128 : *
129 : IF (BNRM2.LT.XPETIT) BNRM2 = ONE
130 : RESID = RNRM2 / BNRM2
131 : IF(IMPR.GT.6) THEN
132 : WRITE(IOIMP,*) 'Résidu initial=',RESID
133 : ENDIF
134 : RESMAX = RESID * RRMAX
135 : LRES.PROG(ITER+1)=RESID
136 : LNMV.LECT(ITER+1)=INMV
137 : IF (RESID.LE.TOL) GOTO 30
138 : C r[tilde]=r(0)
138 : CALL GCOPY(R,RTLD)
139 : * Autre choix
140 : * DO IN=1,RTLD.A(/1)
141 : * RTLD.A(IN) = ONE
142 : * ENDDO
143 : *
144 : 10 CONTINUE
145 : *
146 : * Perform Preconditionned BiConjugate Gradient Stabilized iteration.
147 : *
148 : ITER = ITER + 1
149 : IF(IMPR.GT.7) THEN
150 : WRITE(IOIMP,*) 'ITER=',ITER
151 : ENDIF
152 : *
153 : C rho(i-1)=r[tilde]'r(i-1)
154 : RHO = GDOT(RTLD,R)
155 : * IF ((ITER.LE.10).OR.((ITER/10)*10.EQ.ITER)) THEN
156 : * WRITE(IOIMP,*) 'ITER=',ITER,' RHO=',RHO
157 : * ENDIF
158 : C if rho(i-1)=0 method fails
159 : IF (ABS(RHO).LT.RHOTOL) THEN
160 : RHO=SIGN(XZPREC,RHO)
161 : IWARN=IWARN+1
162 : ENDIF
163 : *
164 : * Compute vector P.
165 : *
166 : IF (ITER.GT.1) THEN
167 : C beta(i-1)=(rho(i-1)/rho(i-2))*
168 : C (alpha(i-1)/omega(i-1))
169 : BETA = ( RHO / RHO1 ) * ( ALPHA / OMEGA )
170 : C p(i)=r(i-1)+beta(i-1)*
171 : C (p(i-1)-omega(i-1)*v(i-1))
171 : CALL GAXPY(-OMEGA,V,P)
171 : CALL GSCAL(BETA,P)
171 : CALL GAXPY(ONE,R,P)
172 : ELSE
173 : C p(1)=r(0)
173 : CALL GCOPY(R,P)
174 : ENDIF
175 : *
176 : * Compute direction adjusting vector PHAT and scalar ALPHA.
177 : *
178 : C
179 : C Preconditionner solve
180 : C Mp[hat]=p(i)
180 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
181 : $ ,PHAT,P)
182 :
183 : C v(i)=Ap[hat]
183 : CALL GMOMV('N',ONE,KMORS,KISA,PHAT,ZERO,V)
184 : INMV=INMV+1
185 : C alpha(i)=rho(i-1)/(r[tilde]'v(i))
186 : TAU = GDOT(RTLD,V)
187 : IF (ABS(TAU).LT.TAUTOL) THEN
188 : TAU=SIGN(XZPREC,TAU)
189 : IWARN=IWARN+1
190 : ENDIF
191 : ALPHA = RHO / TAU
192 : *
193 : * Early check for tolerance.
194 : *
195 : C s=r(i-1)-alpha(i)v(i)
195 : CALL GAXPY(-ALPHA,V,R)
195 : CALL GCOPY(R,S)
196 : RESID = GNRM2(S) / BNRM2
197 : IF(IMPR.GT.6) THEN
198 : WRITE(IOIMP,*) 'Résidu 1 =',RESID
199 : ENDIF
200 : IF (RESID.GT.RESMAX) THEN
201 : IGRAN=IGRAN+1
202 : ELSE
203 : IGRAN=0
204 : ENDIF
205 : IF (RESID.LE.TOL) THEN
206 : C x(i)=x(i-1)+alpha(i)p[hat]
206 : CALL GAXPY(ALPHA,PHAT,X)
207 : LRES.PROG(ITER+1)=RESID
208 : LNMV.LECT(ITER+1)=INMV
209 : GOTO 30
210 : ELSE
211 : *
212 : * Compute stabilizer vector SHAT and scalar OMEGA.
213 : *
214 : C
215 : C Preconditionner solve (ILU(0))
216 : C Ms[hat]=s(i)
216 : CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
217 : $ ,SHAT,S)
218 : C t=As[hat]
218 : CALL GMOMV('N',ONE,KMORS,KISA,SHAT,ZERO,T)
219 : INMV=INMV+1
220 : C
221 : C Il devrait y avoir un test ici ???
222 : C
223 : C omega(i)=t's/t't
224 : OMEGA = GDOT(T,S) / GDOT(T,T)
225 : *
226 : * Compute new solution approximation vector X.
227 : *
228 : C x(i)=x(i-1)+alpha(i)p[hat]
229 : C +omega(i)s[hat]
229 : CALL GAXPY(ALPHA,PHAT,X)
229 : CALL GAXPY(OMEGA,SHAT,X)
230 : *
231 : * Compute residual R, check for tolerance.
232 : *
233 : C r(i)=s-omega(i)t
233 : CALL GAXPY(-OMEGA,T,R)
234 : *
235 : RESID = GNRM2(R) / BNRM2
236 : IF(IMPR.GT.6) THEN
237 : WRITE(IOIMP,*) 'Résidu 2 =',RESID
238 : ENDIF
239 : IF (RESID.GT.RESMAX) THEN
240 : IGRAN=IGRAN+1
241 : ELSE
242 : IGRAN=0
243 : ENDIF
244 : LRES.PROG(ITER+1)=RESID
245 : LNMV.LECT(ITER+1)=INMV
246 : IF (RESID.LE.TOL) GOTO 30
247 : *
248 : * Iteration fails.
249 : *
250 : * IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
251 : IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
252 : IRET = 1
253 : IF (IMPR.GT.0) THEN
254 : C WRITE(IOIMP,*)
255 : C $ 'bcg2.eso : Convergence to tolerance not achieved'
256 : C WRITE(IOIMP,*) 'INMV=',INMV,
257 : C $ ' TOL=',TOL,
258 : C $ ' RESID=',RESID
259 : IF (IWARN.GT.0) THEN
260 : WRITE(IOIMP,*) 'BiCGStab nearly failed ',IWARN
261 : $ ,' times.'
262 : ENDIF
263 : ENDIF
264 : RETURN
265 : ENDIF
266 : *
267 : ENDIF
268 : *
269 : IF (ABS(OMEGA).LT.OMETOL) THEN
270 : OMEGA=SIGN(XZPREC,OMEGA)
271 : IWARN=IWARN+1
272 : ENDIF
273 : *
274 : * Do some more
275 : *
276 : RHO1 = RHO
277 : GOTO 10
278 : *
279 : * Iteration successful
280 : *
281 : 30 CONTINUE
282 : IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
283 : WRITE(IOIMP,*) 'BiCGStab nearly failed ',IWARN
284 : $ ,' times.'
285 : ENDIF
286 : IRET = 0
287 : RETURN
288 : *
289 : * End of BCG2.
290 : *
291 : END
292 :
293 :
294 :
295 :
296 :
297 :
298 :
299 :
300 :
301 :
© Cast3M 2003 - All rights reserved.
Disclaimer