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