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