1 : C BCGSI0 SOURCE GOUNAND 05/06/06 21:15:01 5078
2 : SUBROUTINE BCGSI0(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ JU,JLU,ALU,
5 : $ LRES,
6 : $ R,RTLD,P,PHAT,
7 : $ V,S,SHAT,T,
8 : $ ITER,BRTOL,RESID,IMPR,IRET)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : IMPLICIT INTEGER(I-N)
11 : C***********************************************************************
12 : C NOM : BCGSI0
13 : C DESCRIPTION :
14 : C Résolution d'un système linéaire Ax=b
15 : C par la méthode du bi-gradient conjugué stabilisé
16 : C (Bi-CGSTAB) préconditionné par une factorisation ILU(0)
17 : C (Crout incomplet) ou MILU(0) (Crout inc. modifié relaxé) de A.
18 : C
19 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
20 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
21 : C mél : gounand@semt2.smts.cea.fr
22 : C REFERENCE (bibtex-like) :
23 : C @BOOK{templates,
24 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
25 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
26 : C H. Van der Vorst},
27 : C TITLE={Templates for the Solution of Linear Systems :
28 : C Building Blocks for Iterative Methods},
29 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
30 : C -> URL : http://www.netlib.org/templates/Templates.html
31 : C***********************************************************************
32 : C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines)
33 : C DNRM2, DDOT,DSCAL (functions)
34 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
35 : C PSILU0 : Preconditionner solve
36 : C (Montée-Descente adaptée à ILU(0))
37 : C***********************************************************************
38 : C ENTREES : N, NNZ,
39 : C ROWPTR, COLIND, VAL,
40 : C B, JU, JLU, ALU, BRTOL, IMPR
41 : C ENTREES/SORTIES : X,
42 : C R, RTLD, P, PHAT,
43 : C V, S, SHAT, T,
44 : C ITER, RESID
45 : C SORTIES : LRES, IRET
46 : C CODE RETOUR (IRET) : 0 si ok
47 : C <0 si problème (cf. BRTOL)
48 : C >0 si l'algorithme n'a pas convergé
49 : C N : nombre de degrés de liberté
50 : C NNZ : nombre de valeurs non nulles de la matrice Morse
51 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
52 : C et valeurs de la matrice Morse
53 : C B : vecteur second membre
54 : C ILUVAL : valeurs de la matrice de la factorisation
55 : C ILU(0) ou MILU(0) de A.
56 : C Cette matrice est stockée en Morse, elle a le
57 : C meme profil que A.
58 : C BRTOL : REAL*8
59 : C 'Breakdown tolerance' pour les méthodes de type
60 : C BiCGSTAB. Si un certain produit scalaire de vecteurs
61 : C "direction" est inférieur à cette tolérance, la
62 : C méthode ne peut pas continuer.
63 : C IMPR : niveau d'impression (0..4)
64 : C X : vecteur des inconnues.
65 : C En entrée, contient l'estimation x(0).
66 : C En sortie, la solution trouvée (si la méthode
67 : C a convergé).
68 : C R,RTLD,P,PHAT: vecteurs de travail initialisés et détruits
69 : C V,S,SHAT,T dans prbcgs.
70 : C ITER : type INTEGER.
71 : C En entrée, contient le nombre maximum
72 : C d'itérations à effectuer.
73 : C En sortie, contient le nombre d'itérations
74 : C réellement effectuées.
75 : C RESID : type REAL*8.
76 : C En entrée, la valeur maximale du résidu à
77 : C convergence de l'algorithme mesurée comme suit :
78 : C norme[L2](b - A*x) / norme[L2]( b )
79 : C En sortie, la valeur effective de cette mesure.
80 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
81 : C contient l'historique de la valeur de RESID en
82 : C fonction du nombre d'itérations effectuées.
83 : C
84 : C***********************************************************************
85 : C VERSION : 20/12/99
86 : C HISTORIQUE : v1, 02/04/98, création
87 : C HISTORIQUE : 20/12/99 : Interfaçage avec le nouveau psilu0
88 : C (factorisations incomplètes stockées au format MSR (Modified
89 : C Sparse Row).
90 : C HISTORIQUE :
91 : C***********************************************************************
92 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
93 : C en cas de modification de ce sous-programme afin de faciliter
94 : C la maintenance !
95 : C***********************************************************************
96 : *
97 : * .. Includes .. a commenter si mise au point ok
98 : -INC CCREEL
99 : -INC CCOPTIO
100 : INTEGER IMPR
101 : * .. Scalar Arguments ..
102 : INTEGER N,NNZ
103 : INTEGER ITER,IRET
104 : REAL*8 BRTOL,RESID
105 : * ..
106 : * .. Array Arguments ..
107 : * .. Matrices stockées en Morse
108 : INTEGER ROWPTR(N+1)
109 : INTEGER COLIND(NNZ)
110 : REAL*8 VAL(NNZ)
111 : INTEGER JU(N)
112 : INTEGER JLU(NNZ+1)
113 : REAL*8 ALU(NNZ)
114 : * .. Vecteurs
115 : REAL*8 B(N)
116 : REAL*8 X(N)
117 : REAL*8 LRES(ITER)
118 : REAL*8 R(N),RTLD(N),P(N),PHAT(N)
119 : REAL*8 V(N),S(N),SHAT(N),T(N)
120 : *
121 : *
122 : * .. Variables locales
123 : * .. Parameters
124 : REAL*8 ZERO,ONE
125 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
126 : * ..
127 : INTEGER MAXIT,III
128 : REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
129 : REAL*8 RHOTOL,OMETOL
130 : REAL*8 DDOT,DNRM2
131 : * ..
132 : * .. External subroutines and functions..
133 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL
134 : * ..
135 : * .. Executable Statements ..
136 : *
137 : *
138 : IRET = 0
139 : IWARN = 0
140 : MAXIT = ITER
141 : TOL = RESID
142 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgilu0.eso'
143 : *
144 : * Set parameter tolerances.
145 : *
146 : RHOTOL = BRTOL
147 : OMETOL = BRTOL
148 : *
149 : * Set initial residual.
150 : *
151 : C r(0)=b
151 : CALL DCOPY(N,B,1,R,1)
152 : IF (DNRM2(N,X,1).NE.ZERO) THEN
153 : C r(0)=b-Ax
153 : CALL DMOMV ('No transpose',N,-ONE,
154 : $ NNZ,ROWPTR,COLIND,VAL,
155 : $ X,ONE,R)
156 : ENDIF
157 : *
158 : ITER = 0
159 :
160 : BNRM2 = DNRM2(N,B,1)
161 : IF(IMPR.GT.6) THEN
162 : WRITE(IOIMP,*) '||B||=',BNRM2
163 : ENDIF
164 :
165 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
166 : RESID = DNRM2(N,R,1) / BNRM2
167 : * RESID = DNRM2(N,R,1)
168 : IF(IMPR.GT.6) THEN
169 : WRITE(IOIMP,*) 'Résidu initial=',RESID
170 : ENDIF
171 :
172 : IF (RESID.LE.TOL) GOTO 30
173 : C r[tilde]=r(0)
173 : *!!! CALL DCOPY(N,R,1,RTLD,1)
174 : * Autre choix
175 : DO IN=1,N
176 : RTLD(IN) = ONE
177 : ENDDO
178 : *
179 : 10 CONTINUE
180 : *
181 : * Perform BiConjugate Gradient Stabilized iteration.
182 : *
183 : ITER = ITER + 1
184 : IF(IMPR.GT.7) THEN
185 : WRITE(IOIMP,*) 'ITER=',ITER
186 : ENDIF
187 : *
188 : C rho(i-1)=r[tilde]'r(i-1)
189 : RHO = DDOT(N,RTLD,1,R,1)
190 : C if rho(i-1)=0 method fails
191 : * IF (ABS(RHO).LT.RHOTOL) GOTO 25
192 : IF (ABS(RHO).LT.RHOTOL) THEN
193 : RHO=XZPREC
194 : IWARN=IWARN+1
195 : ENDIF
196 : *
197 : * Compute vector P.
198 : *
199 : IF (ITER.GT.1) THEN
200 : C beta(i-1)=(rho(i-1)/rho(i-2))*
201 : C (alpha(i-1)/omega(i-1))
202 : BETA = ( RHO / RHO1 ) * ( ALPHA / OMEGA )
203 : C p(i)=r(i-1)+beta(i-1)*
204 : C (p(i-1)-omega(i-1)*v(i-1))
204 : CALL DAXPY(N,-OMEGA,V,1,P,1)
204 : CALL DSCAL(N,BETA,P,1)
204 : CALL DAXPY(N,ONE,R,1,P,1)
205 : ELSE
206 : C p(1)=r(0)
206 : CALL DCOPY(N,R,1,P,1)
207 : ENDIF
208 : *
209 : * Compute direction adjusting vector PHAT and scalar ALPHA.
210 : *
211 : C
212 : C Preconditionner solve (ILU(0))
213 : C Mp[hat]=p(i)
213 : CALL PSILU0(N,NNZ,JU,JLU,ALU,
214 : $ PHAT,P)
215 : C v(i)=Ap[hat]
215 : CALL DMOMV ('No transpose',N,ONE,
216 : $ NNZ,ROWPTR,COLIND,VAL,
217 : $ PHAT,ZERO,V)
218 : C alpha(i)=rho(i-1)/(r[tilde]'v(i))
219 : ALPHA = RHO / DDOT(N,RTLD,1,V,1)
220 : *
221 : * Early check for tolerance.
222 : *
223 : C s=r(i-1)-alpha(i)v(i)
223 : CALL DAXPY(N,-ALPHA,V,1,R,1)
223 : CALL DCOPY(N,R,1,S,1)
224 : RESID = DNRM2(N,S,1) / BNRM2
225 : IF(IMPR.GT.6) THEN
226 : WRITE(IOIMP,*) 'Résidu 1 =',RESID
227 : ENDIF
228 : IF (RESID.LE.TOL) THEN
229 : C x(i)=x(i-1)+alpha(i)p[hat]
229 : CALL DAXPY(N,ALPHA,PHAT,1,X,1)
230 : LRES(ITER)=RESID
231 : GOTO 30
232 : ELSE
233 :
234 : *
235 : * Compute stabilizer vector SHAT and scalar OMEGA.
236 : *
237 : C
238 : C Preconditionner solve (ILU(0))
239 : C Ms[hat]=s(i)
239 : CALL PSILU0(N,NNZ,JU,JLU,ALU,
240 : $ SHAT,S)
241 : C t=As[hat]
241 : CALL DMOMV ('No transpose',N,ONE,
242 : $ NNZ,ROWPTR,COLIND,VAL,
243 : $ SHAT,ZERO,T)
244 : C omega(i)=t's/t't
245 : OMEGA = DDOT(N,T,1,S,1) /
246 : $ DDOT(N,T,1,T,1)
247 : *
248 : * Compute new solution approximation vector X.
249 : *
250 : C x(i)=x(i-1)+alpha(i)p[hat]
251 : C +omega(i)s[hat]
251 : CALL DAXPY(N,ALPHA,PHAT,1,X,1)
251 : CALL DAXPY(N,OMEGA,SHAT,1,X,1)
252 : IF(IMPR.GT.7) THEN
253 : WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
254 : ENDIF
255 :
256 : *
257 : * Compute residual R, check for tolerance.
258 : *
259 : C r(i)=s-omega(i)t
259 : CALL DAXPY(N,-OMEGA,T,1,R,1)
260 : *
261 : RESID = DNRM2(N,R,1) / BNRM2
262 : IF(IMPR.GT.6) THEN
263 : WRITE(IOIMP,*) 'Résidu 2 =',RESID
264 : ENDIF
265 : LRES(ITER)=RESID
266 : IF (RESID.LE.TOL) GOTO 30
267 : *
268 : * Iteration fails.
269 : *
270 : IF (ITER.EQ.MAXIT) THEN
271 : IRET = 1
272 : WRITE(IOIMP,*)
273 : $ 'bcgsi0.eso : Convergence to tolerance not achieved'
274 : WRITE(IOIMP,*) 'ITER=',ITER,
275 : $ ' TOL=',TOL,
276 : $ ' RESID=',RESID
277 : RETURN
278 : ENDIF
279 : *
280 : ENDIF
281 : *
282 : * IF (ABS(OMEGA).LT.OMETOL) GOTO 25
283 : IF (ABS(OMEGA).LT.OMETOL) THEN
284 : OMEGA=XZPREC
285 : IWARN=IWARN+1
286 : ENDIF
287 : *
288 : * Do some more
289 : *
290 : RHO1 = RHO
291 : GOTO 10
292 : *
293 : * A breakdown in the method has occured.
294 : *
295 : 25 CONTINUE
296 : IF (ABS(RHO).LT.RHOTOL) THEN
297 : IRET=-1
298 : ELSEIF(ABS(OMEGA).LT.OMETOL) THEN
299 : IRET=-2
300 : ENDIF
301 : RETURN
302 : *
303 : * Iteration successful
304 : *
305 : 30 CONTINUE
306 : IF (IWARN.GT.0) THEN
307 : WRITE(IOIMP,*) 'BiCGStab a failli planter ',IWARN,' fois.'
308 : ENDIF
309 : IRET = 0
310 : RETURN
311 : *
312 : * End of BCGSI0.
313 : *
314 : END
315 :
316 :
317 :
318 :
319 :
320 :
321 :
© Cast3M 2003 - All rights reserved.
Disclaimer