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