1 : C BCGSDI SOURCE CHAT 05/01/12 21:37:06 5004
2 : SUBROUTINE BCGSDI(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ DINV,
5 : $ LRES,
6 : $ R,RTLD,P,PHAT,
7 : $ V,S,SHAT,T,
8 : $ ITER,BRTOL,RESID,IMPR,IRET)
9 : C***********************************************************************
10 : C NOM : BCGSDI
11 : C DESCRIPTION :
12 : C Résolution d'un système linéaire Ax=b
13 : C par la méthode du bi-gradient conjugué stabilisé
14 : C (Bi-CGSTAB) préconditionné par la factorisation D-ILU de A.
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***********************************************************************
29 : C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines)
30 : C DNRM2, DDOT,DSCAL (functions)
31 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
32 : C PSDILU : Preconditionner solve
33 : C (Montée-Descente adaptée à D-ILU)
34 : C***********************************************************************
35 : C ENTREES : N, NNZ,
36 : C ROWPTR, COLIND, VAL,
37 : C B, DINV, BRTOL, IMPR
38 : C ENTREES/SORTIES : X,
39 : C R, RTLD, P, PHAT,
40 : C V, S, SHAT, T,
41 : C ITER, RESID
42 : C SORTIES : LRES, IRET
43 : C CODE RETOUR (IRET) : 0 si ok
44 : C <0 si problème (cf. BRTOL)
45 : C >0 si l'algorithme n'a pas convergé
46 : C N : nombre de degrés de liberté
47 : C NNZ : nombre de valeurs non nulles de la matrice Morse
48 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
49 : C et valeurs de la matrice Morse
50 : C B : vecteur second membre
51 : C DINV : vecteur contenant les inverses des pivots
52 : C de la factorisation D-ILU de A.
53 : C BRTOL : REAL*8
54 : C 'Breakdown tolerance' pour les méthodes de type
55 : C BiCGSTAB. Si un certain produit scalaire de vecteurs
56 : C "direction" est inférieur à cette tolérance, la
57 : C méthode ne peut pas continuer.
58 : C IMPR : niveau d'impression (0..4)
59 : C X : vecteur des inconnues.
60 : C En entrée, contient l'estimation x(0).
61 : C En sortie, la solution trouvée (si la méthode
62 : C a convergé).
63 : C R,RTLD,P,PHAT: vecteurs de travail initialisés et détruits
64 : C V,S,SHAT,T dans prbcgs.
65 : C ITER : type INTEGER.
66 : C En entrée, contient le nombre maximum
67 : C d'itérations à effectuer.
68 : C En sortie, contient le nombre d'itérations
69 : C réellement effectuées.
70 : C RESID : type REAL*8.
71 : C En entrée, la valeur maximale du résidu à
72 : C convergence de l'algorithme mesurée comme suit :
73 : C norme[L2](b - A*x) / norme[L2]( b )
74 : C En sortie, la valeur effective de cette mesure.
75 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
76 : C contient l'historique de la valeur de RESID en
77 : C fonction du nombre d'itérations effectuées.
78 : C
79 : C***********************************************************************
80 : C VERSION : v1, 02/04/98, version initiale
81 : C HISTORIQUE : v1, 02/04/98, création
82 : C HISTORIQUE :
83 : C HISTORIQUE :
84 : C***********************************************************************
85 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
86 : C en cas de modification de ce sous-programme afin de faciliter
87 : C la maintenance !
88 : C***********************************************************************
89 : *
90 : * .. Includes .. a commenter si mise au point ok
91 : IMPLICIT INTEGER(I-N)
92 : -INC CCOPTIO
93 : INTEGER IMPR
94 : * .. Scalar Arguments ..
95 : INTEGER N,NNZ
96 : INTEGER ITER,IRET
97 : REAL*8 BRTOL,RESID
98 : * ..
99 : * .. Array Arguments ..
100 : * .. Matrices stockées en Morse
101 : INTEGER ROWPTR(N+1)
102 : INTEGER COLIND(NNZ)
103 : REAL*8 VAL(NNZ)
104 : * .. Vecteurs
105 : REAL*8 B(N)
106 : REAL*8 X(N)
107 : REAL*8 DINV(N)
108 : REAL*8 LRES(ITER)
109 : REAL*8 R(N),RTLD(N),P(N),PHAT(N)
110 : REAL*8 V(N),S(N),SHAT(N),T(N)
111 : *
112 : *
113 : * .. Variables locales
114 : * .. Parameters
115 : REAL*8 ZERO,ONE
116 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
117 : * ..
118 : INTEGER MAXIT
119 : REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
120 : REAL*8 RHOTOL,OMETOL
121 : REAL*8 DDOT,DNRM2
122 : * ..
123 : * .. External subroutines and functions..
124 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL
125 : * ..
126 : * .. Executable Statements ..
127 : *
128 : *
129 : IRET = 0
130 : MAXIT = ITER
131 : TOL = RESID
132 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgsdi.eso'
133 : *
134 : * Set parameter tolerances.
135 : *
136 : RHOTOL = BRTOL
137 : OMETOL = BRTOL
138 : *
139 : * Set initial residual.
140 : *
141 : C r(0)=b
141 : CALL DCOPY(N,B,1,R,1)
142 : IF (DNRM2(N,X,1).NE.ZERO) THEN
143 : C r(0)=b-Ax
143 : CALL DMOMV ('No transpose',N,-ONE,
144 : $ NNZ,ROWPTR,COLIND,VAL,
145 : $ X,ONE,R)
146 : ENDIF
147 : *
148 : ITER = 0
149 :
150 : BNRM2 = DNRM2(N,B,1)
151 : IF(IMPR.GT.5) THEN
152 : WRITE(IOIMP,*) '||B||=',BNRM2
153 : ENDIF
154 :
155 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
156 : RESID = DNRM2(N,R,1) / BNRM2
157 : * RESID = DNRM2(N,R,1)
158 : IF(IMPR.GT.5) THEN
159 : WRITE(IOIMP,*) 'Résidu initial=',RESID
160 : ENDIF
161 :
162 : IF (RESID.LE.TOL) GOTO 30
163 : C r[tilde]=r(0)
163 : CALL DCOPY(N,R,1,RTLD,1)
164 : *
165 : 10 CONTINUE
166 : *
167 : * Perform BiConjugate Gradient Stabilized iteration.
168 : *
169 : ITER = ITER + 1
170 : IF(IMPR.GT.6) THEN
171 : WRITE(IOIMP,*) 'ITER=',ITER
172 : ENDIF
173 : *
174 : C rho(i-1)=r[tilde]'r(i-1)
175 : RHO = DDOT(N,RTLD,1,R,1)
176 : IF(IMPR.GT.6) THEN
177 : WRITE(IOIMP,*) 'RHO=',RHO
178 : ENDIF
179 : C if rho(i-1)=0 method fails
180 : IF (ABS(RHO).LT.RHOTOL) GOTO 25
181 : *
182 : * Compute vector P.
183 : *
184 : IF (ITER.GT.1) THEN
185 : C beta(i-1)=(rho(i-1)/rho(i-2))*
186 : C (alpha(i-1)/omega(i-1))
187 : BETA = ( RHO / RHO1 ) * ( ALPHA / OMEGA )
188 : C p(i)=r(i-1)+beta(i-1)*
189 : C (p(i-1)-omega(i-1)*v(i-1))
189 : CALL DAXPY(N,-OMEGA,V,1,P,1)
189 : CALL DSCAL(N,BETA,P,1)
189 : CALL DAXPY(N,ONE,R,1,P,1)
190 : ELSE
191 : C p(1)=r(0)
191 : CALL DCOPY(N,R,1,P,1)
192 : ENDIF
193 : *
194 : * Compute direction adjusting vector PHAT and scalar ALPHA.
195 : *
196 : C
197 : C Preconditionner solve (D-ILU)
198 : C
199 : C Mp[hat]=p(i)
199 : CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
200 : $ DINV,
201 : $ PHAT,P)
202 : C v(i)=Ap[hat]
202 : CALL DMOMV ('No transpose',N,ONE,
203 : $ NNZ,ROWPTR,COLIND,VAL,
204 : $ PHAT,ZERO,V)
205 : C alpha(i)=rho(i-1)/(r[tilde]'v(i))
206 : ALPHA = RHO / DDOT(N,RTLD,1,V,1)
207 : *
208 : * Early check for tolerance.
209 : *
210 : C s=r(i-1)-alpha(i)v(i)
210 : CALL DAXPY(N,-ALPHA,V,1,R,1)
210 : CALL DCOPY(N,R,1,S,1)
211 : RESID = DNRM2(N,S,1) / BNRM2
212 : IF(IMPR.GT.6) THEN
213 : WRITE(IOIMP,*) 'Résidu 1 =',RESID
214 : ENDIF
215 : IF (RESID.LE.TOL) THEN
216 : C x(i)=x(i-1)+alpha(i)p[hat]
216 : CALL DAXPY(N,ALPHA,PHAT,1,X,1)
217 : LRES(ITER)=RESID
218 : GOTO 30
219 : ELSE
220 :
221 : *
222 : * Compute stabilizer vector SHAT and scalar OMEGA.
223 : *
224 : C
225 : C Preconditionner solve (D-ILU)
226 : C
227 : C Ms[hat]=s(i)
227 : CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
228 : $ DINV,
229 : $ SHAT,S)
230 :
231 : C t=As[hat]
231 : CALL DMOMV ('No transpose',N,ONE,
232 : $ NNZ,ROWPTR,COLIND,VAL,
233 : $ SHAT,ZERO,T)
234 : C omega(i)=t's/t't
235 : OMEGA = DDOT(N,T,1,S,1) /
236 : $ DDOT(N,T,1,T,1)
237 : IF(IMPR.GT.6) THEN
238 : WRITE(IOIMP,*) 'OMEGA=',OMEGA
239 : ENDIF
240 : *
241 : * Compute new solution approximation vector X.
242 : *
243 : C x(i)=x(i-1)+alpha(i)p[hat]
244 : C +omega(i)s[hat]
244 : CALL DAXPY(N,ALPHA,PHAT,1,X,1)
244 : CALL DAXPY(N,OMEGA,SHAT,1,X,1)
245 : IF(IMPR.GT.7) THEN
246 : WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
247 : ENDIF
248 :
249 : *
250 : * Compute residual R, check for tolerance.
251 : *
252 : C r(i)=s-omega(i)t
252 : CALL DAXPY(N,-OMEGA,T,1,R,1)
253 : *
254 : RESID = DNRM2(N,R,1) / BNRM2
255 : IF(IMPR.GT.6) THEN
256 : WRITE(IOIMP,*) 'Résidu 2 =',RESID
257 : ENDIF
258 : LRES(ITER)=RESID
259 : IF (RESID.LE.TOL) GOTO 30
260 : *
261 : * Iteration fails.
262 : *
263 : IF (ITER.EQ.MAXIT) THEN
264 : IRET = 1
265 : WRITE(IOIMP,*)
266 : $ 'bcgsdi.eso : Convergence to tolerance not achieved'
267 : WRITE(IOIMP,*) 'ITER=',ITER
268 : WRITE(IOIMP,*) 'TOL=',TOL
269 : WRITE(IOIMP,*) 'RESID=',RESID
270 : RETURN
271 : ENDIF
272 : *
273 : ENDIF
274 : *
275 : IF (ABS(OMEGA).LT.OMETOL) GOTO 25
276 : *
277 : * Do some more
278 : *
279 : RHO1 = RHO
280 : GOTO 10
281 : *
282 : * A breakdown in the method has occured.
283 : *
284 : 25 CONTINUE
285 : IF (ABS(RHO).LT.RHOTOL) THEN
286 : IRET=-1
287 : ELSEIF(ABS(OMEGA).LT.OMETOL) THEN
288 : IRET=-2
289 : ENDIF
290 : RETURN
291 : *
292 : * Iteration successful
293 : *
294 : 30 CONTINUE
295 : IRET = 0
296 : RETURN
297 : *
298 : * End of BCGSDI.
299 : *
300 : END
301 :
302 :
303 :
304 :
© Cast3M 2003 - All rights reserved.
Disclaimer