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