1 : C PRBCGS SOURCE GOUNAND 11/07/21 21:15:38 7046
2 : SUBROUTINE PRBCGS(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,INCX,ITER,
3 : $ RESID,KPREC,
4 : $ BRTOL,IMPR,IRET)
5 : IMPLICIT REAL*8 (A-H,O-Z)
6 : IMPLICIT INTEGER (I-N)
7 : C***********************************************************************
8 : C NOM : PRBCGS
9 : C DESCRIPTION :
10 : C Préparation de la résolution d'un système linéaire Ax=b
11 : C par une méthode de bi-gradient conjugué stabilisé
12 : C (Bi-CGSTAB <~> Bi-CG + GMRES(1))
13 : C préconditionnée ou non.
14 : C
15 : C A peut etre symétrique (mais mieux vaut utiliser
16 : C le gradient conjugué)
17 : C ou non symétrique (mais la méthode peut ne pas converger
18 : C voir paramètre BRTOL en particulier).
19 : C
20 : C Les 4 préconditionnements disponibles pour cette méthode
21 : C sont :
22 : C * Jacobi (diagonal)
23 : C * D-ILU (Diagonal Incomplete LU factorization)
24 : C * ILU(0) (Incomplete LU factorization
25 : C of level zero ie Choleski incomplet)
26 : C * MILU(0) (Modified ILU(0)) relaxé
27 : C 5 : préconditionnement ILUT (dual truncation)
28 : C 6 : préconditionnement ILUT2 (une variante du
29 : C précédent qui remplit mieux la mémoire et
30 : C fonctionne mieux quelquefois)
31 : C
32 : C Ces préconditionneurs sont supposés déjà construits (ie
33 : C factorisés), voir les subroutines MEIDIA, PRDILU, PRILU0,
34 : C PRMILU et PRILUT.
35 : C
36 : C Ce sous-programme est en fait une interface aux sous-programmes :
37 : C bcgs, bcgsd, bcgsdi, bcgsi0
38 : C qui sont en Fortran presque pur (pour raison de rapidité)
39 : C et effectuent la résolution du système linéaire
40 : C proprement dite.
41 : C
42 : C LANGAGE : ESOPE
43 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
44 : C mél : gounand@semt2.smts.cea.fr
45 : C REFERENCE (bibtex-like) :
46 : C @BOOK{templates,
47 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
48 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
49 : C H. Van der Vorst},
50 : C TITLE={Templates for the Solution of Linear Systems :
51 : C Building Blocks for Iterative Methods},
52 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
53 : C -> URL : http://www.netlib.org/templates/Templates.html
54 : C***********************************************************************
55 : C APPELES (Calcul) : BCGS, BCGSD, BCGSDI, BCGSI0
56 : C***********************************************************************
57 : C ENTREES : MATRIK, MAPREC, KPREC, BRTOL, IMPR
58 : C ENTREES/SORTIES : INCX, ITER, RESID
59 : C SORTIES : LRES, IRET
60 : C CODE RETOUR (IRET) : 0 si ok
61 : C <0 si problème
62 : C >0 si l'algorithme appelé n'a pas convergé
63 : C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
64 : C on pioche dedans les informations nécessaires
65 : C (différents pointeurs, nb. de ddl...)
66 : C MAPREC : pointeur sur segment MATRIK de l'include SMMATRIK
67 : C on va l'utiliser comme préconditionneur.
68 : C KPREC : type de préconditionnement à effectuer
69 : C 0 : pas de préconditionnement
70 : C 1 : préconditionnement Jacobi (diagonal)
71 : C 2 : préconditionnement D-ILU
72 : C 3 : préconditionnement ILU(0) (Choleski incomplet)
73 : C 4 : préconditionnement MILU(0)
74 : C (Choleski incomplet modifié)
75 : C 5 : préconditionnement ILUT (dual truncation)
76 : C 6 : préconditionnement ILUT2 (une variante du
77 : C précédent qui remplit mieux la mémoire et
78 : C fonctionne mieux quelquefois)
79 : C BRTOL : REAL*8
80 : C 'Breakdown tolerance' pour les méthodes de type
81 : C BiCGSTAB. Si un certain produit scalaire de vecteurs
82 : C "direction" est inférieur à cette tolérance, la
83 : C méthode ne peut pas continuer.
84 : C IMPR : niveau d'impression (0..4)
85 : C INCX : pointeur sur segment IZA de l'include SMMATRIK
86 : C C'est le vecteur des inconnues.
87 : C En entrée, contient l'estimation x(0).
88 : C En sortie, la solution trouvée (si la méthode
89 : C a convergé).
90 : C ITER : type INTEGER.
91 : C En entrée, contient le nombre maximum
92 : C d'itérations à effectuer.
93 : C En sortie, contient le nombre d'itérations
94 : C réellement effectuées.
95 : C RESID : type REAL*8.
96 : C En entrée, la valeur maximale du résidu à
97 : C convergence de l'algorithme mesurée comme suit :
98 : C norme[L2](b - A*x) / norme[L2]( b )
99 : C En sortie, la valeur effective de cette mesure.
100 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
101 : C contient l'historique de la valeur de RESID en
102 : C fonction du nombre d'itérations effectuées.
103 : C
104 : C***********************************************************************
105 : C VERSION : v1, 02/04/98, version initiale
106 : C HISTORIQUE : v1, 02/04/98, création
107 : C HISTORIQUE : 09/02/99 : ajout préconditionneur <> MATRIK
108 : C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau bcgsi0
109 : C (factorisations incomplètes stockées au format MSR (Modified Sparse
110 : C Row) (voir la doc de Sparskit version 2+)
111 : C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
112 : C HISTORIQUE : 22/03/00 : ajout des préconditionneurs ILUT
113 : C HISTORIQUE : 09/04/04 : ajout des préconditionneurs ILUTP
114 : C HISTORIQUE :
115 : C***********************************************************************
116 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
117 : C en cas de modification de ce sous-programme afin de faciliter
118 : C la maintenance !
119 : C***********************************************************************
120 : *
121 : * .. Includes et pointeurs associés ..
122 : -INC CCOPTIO
123 : -INC SMLREEL
124 : INTEGER JG
125 : POINTEUR LRES.MLREEL
126 : -INC SMMATRIK
127 : POINTEUR MAPREC.MATRIK
128 : POINTEUR KMORS.PMORS
129 : POINTEUR KISA.IZA
130 : POINTEUR KS2B.IZA
131 : POINTEUR INCX.IZA
132 : POINTEUR INVDIA.IZA
133 : POINTEUR ILUM.PMORS
134 : POINTEUR ILUI.IZA
135 : * .. Parameters
136 : * This is a breakdown tolerance for BiCGSTAB-type method
137 : REAL*8 BRTOL
138 : * .. Work Vectors for BCGS..
139 : POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
140 : POINTEUR IV.IZA,IS.IZA,ISHAT.IZA,IT.IZA
141 : *STAT-INC SMSTAT
142 : * .. Scalar Arguments ..
143 : INTEGER ITER, KPREC, IMPR, IRET
144 : REAL*8 RESID
145 : * ..
146 : * Vars reqd for STOPTEST2
147 : * REAL*8 TOL, BNRM2
148 : * ..
149 : * .. External subroutines ..
150 : * EXTERNAL STOPTEST2
151 : INTEGER NBVA,NJA,NTT,NTTPRE
152 : * ..
153 : * .. Executable Statements ..
154 : *
155 : IRET = 0
156 : *
157 : * On récupère les paramètres
158 : *
159 :
160 : SEGACT MATRIK
161 : SEGACT MAPREC
162 : IF (KSYM.EQ.0) THEN
163 : IF (IMPR.GT.0) THEN
164 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
165 : $ 'Use a Conjugate Gradient instead !'
166 : C IRET=-2
167 : C GOTO 9999
168 : ENDIF
169 : ENDIF
170 : * Pour le préconditionneur
171 : ILUM =MAPREC.KIDMAT(6)
172 : ILUI =MAPREC.KIDMAT(7)
173 : IDMAT=MAPREC.KIDMAT(1)
174 : SEGACT IDMAT
175 : INVDIA=IDIAG
176 : SEGDES IDMAT
177 : SEGACT KMORS
178 : NTT =KMORS.IA(/1)-1
179 : NJA =KMORS.JA(/1)
180 : SEGACT KISA
181 : SEGACT KS2B
182 : SEGACT INCX*MOD
183 : C Paramètres des préconditionnements diagonaux et D-ILU
184 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
185 : C Est-il compatible avec la matrice ?
186 : SEGACT INVDIA
187 : NTTPRE=INVDIA.A(/1)
188 : IF (NTTPRE.NE.NTT) THEN
189 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
190 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
191 : WRITE(IOIMP,*) 'NTT=',NTT
192 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
193 : IRET=-2
194 : GOTO 9999
195 : ENDIF
196 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2
197 : C ,ilutp, ilutpg, ilutpg2
198 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
199 : SEGACT ILUM
200 : SEGACT ILUI
201 : * SEGPRT,ILUM
202 : * SEGPRT,ILUI
203 : NTTPRE=ILUM.IA(/1)
204 : IF (NTTPRE.NE.NTT) THEN
205 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
206 : $ 'ne sont pas compatibles...'
207 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
208 : IRET=-2
209 : GOTO 9999
210 : ENDIF
211 : ELSE
212 : INVDIA=0
213 : ENDIF
214 : C
215 : C Initialisation de l'historique de convergence
216 : C
217 : JG=ITER
218 : SEGINI LRES
219 : C
220 : C
221 : C Initialisation des vecteurs de travail pour BiCGSTAB
222 : C
222 : *STAT CALL INMSTA(MSTAT,0)
223 : NBVA=NTT
224 : SEGINI IR,IRTLD,IP,IPHAT
225 : SEGINI IV,IS,ISHAT,IT
225 : *STAT CALL PRMSTA(' PRBCGS : objets temporaires',MSTAT,1)
225 : *STAT CALL SUMSTA(MSTAT,0)
226 : C
227 : IF (KPREC.EQ.0) THEN
228 : C
229 : C Le BiCGSTAB non préconditionné (fortran pur)
230 : C
230 : CALL BCGS(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
231 : $ KS2B.A,INCX.A,
232 : $ LRES.PROG,
233 : $ IR.A,IRTLD.A,IP.A,IPHAT.A,
234 : $ IV.A,IS.A,ISHAT.A,IT.A,
235 : $ ITER,BRTOL,RESID,IMPR,IRET)
236 : ELSE IF (KPREC.EQ.1) THEN
237 : C
238 : C Le gradient conjugué préconditionné par la diagonale (fortran pur)
239 : C
239 : CALL BCGSD(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
240 : $ KS2B.A,INCX.A,
241 : $ INVDIA.A,
242 : $ LRES.PROG,
243 : $ IR.A,IRTLD.A,IP.A,IPHAT.A,
244 : $ IV.A,IS.A,ISHAT.A,IT.A,
245 : $ ITER,BRTOL,RESID,IMPR,IRET)
246 : ELSE IF (KPREC.EQ.2) THEN
247 : C
248 : C Le gradient conjugué préconditionné par une factorisation
249 : C incomplete D-ILU (fortran pur)
250 : C
250 : CALL BCGSDI(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
251 : $ KS2B.A,INCX.A,
252 : $ INVDIA.A,
253 : $ LRES.PROG,
254 : $ IR.A,IRTLD.A,IP.A,IPHAT.A,
255 : $ IV.A,IS.A,ISHAT.A,IT.A,
256 : $ ITER,BRTOL,RESID,IMPR,IRET)
257 : ELSE IF (KPREC.GE.3.AND.KPREC.LE.9) THEN
258 : C
259 : C Le gradient conjugué préconditionné par une factorisation
260 : C incomplete (fortran pur)
261 : C
261 : CALL BCGSI0(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
262 : $ KS2B.A,INCX.A,
263 : $ ILUM.IA,ILUM.JA,ILUI.A,
264 : $ LRES.PROG,
265 : $ IR.A,IRTLD.A,IP.A,IPHAT.A,
266 : $ IV.A,IS.A,ISHAT.A,IT.A,
267 : $ ITER,BRTOL,RESID,IMPR,IRET)
268 : ELSE
269 : WRITE (IOIMP,*) 'KPREC out of range (0..9)'
270 : IRET=-4
271 : GOTO 9999
272 : ENDIF
273 : C
274 : C Désactivation
275 : C
276 : SEGSUP IR,IRTLD,IP,IPHAT
277 : SEGSUP IV,IS,ISHAT,IT
278 : JG=ITER
279 : SEGADJ LRES
280 : SEGDES LRES
281 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
282 : SEGDES INVDIA
283 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
284 : SEGDES ILUI
285 : SEGDES ILUM
286 : ENDIF
287 : SEGDES INCX
288 : SEGDES KS2B
289 : SEGDES KISA
290 : SEGDES KMORS
291 : SEGDES MAPREC
292 : SEGDES MATRIK
293 : C
294 : C A breakdown has occured in the BiCGSTAB method
295 : C
296 : IF (IRET.LT.0) GOTO 9999
297 : *
298 : * Normal termination
299 : *
300 : RETURN
301 : *
302 : * Format handling
303 : *
304 : 1002 FORMAT(10(1X,1PE11.4))
305 : *
306 : * Error handling
307 : *
308 : 9999 CONTINUE
309 : WRITE(IOIMP,*) 'An error was detected in prbcgs.eso'
310 : RETURN
311 : *
312 : * End of PRBCGS
313 : *
314 : END
315 :
316 :
317 :
318 :
319 :
320 :
321 :
322 :
323 :
324 :
© Cast3M 2003 - All rights reserved.
Disclaimer