1 : C PRGMR SOURCE GOUNAND 11/07/21 21:15:43 7046
2 : SUBROUTINE PRGMR(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,INCX,
3 : $ ITER,RESID,KPREC,
4 : $ RESTRT,IMPR,IRET)
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : C***********************************************************************
8 : C NOM : PRGMR
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 GMRES(m) (restarted Generalized Minimal
12 : C Residual) préconditionnée ou non.
13 : C
14 : C A peut etre symétrique (mais mieux vaut utiliser
15 : C le gradient conjugué)
16 : C ou non symétrique (le résidu b-Ax est toujours DECROISSANT
17 : C mais peut éventuellement stagner si
18 : C m=RESTRT est trop faible).
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 gmr, gmrd, gmrdi, gmri0
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) : GMR, GMRD, GMRDI, GMRI0
56 : C***********************************************************************
57 : C ENTREES : MATRIK, MAPREC, KPREC, RESTRT, 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) (Crout incomplet)
73 : C 4 : préconditionnement MILU(0)
74 : C (Crout 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 RESTRT : INTEGER
80 : C paramètre de redémarrage (restart) (i.e m) pour
81 : C la méthode GMRES(m).
82 : C En gros, on 'minimise' le résidu par rapport
83 : C à ITER<(m=RESTRT) vecteurs orthogonaux 2 à 2.
84 : C A ITER=multiple(m), on réinitialise le tableau
85 : C de vecteurs orthogonaux à zéro.
86 : C Plus m est grand, plus la convergence est bonne
87 : C et plus l'occupation mémoire et le temps de calcul
88 : C moyen pour une itération sont élevés.
89 : C A contrario, si m est trop petit, on peut avoir
90 : C stagnation du résidu.
91 : C IMPR : niveau d'impression (0..4)
92 : C INCX : pointeur sur segment IZA de l'include SMMATRIK
93 : C C'est le vecteur des inconnues.
94 : C En entrée, contient l'estimation x(0).
95 : C En sortie, la solution trouvée (si la méthode
96 : C a convergé).
97 : C ITER : type INTEGER.
98 : C En entrée, contient le nombre maximum
99 : C d'itérations à effectuer.
100 : C En sortie, contient le nombre d'itérations
101 : C réellement effectuées.
102 : C RESID : type REAL*8.
103 : C En entrée, la valeur maximale du résidu à
104 : C convergence de l'algorithme mesurée comme suit :
105 : C norme[L2](b - A*x) / norme[L2]( b )
106 : C En sortie, la valeur effective de cette mesure.
107 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
108 : C contient l'historique de la valeur de RESID en
109 : C fonction du nombre d'itérations effectuées.
110 : C
111 : C***********************************************************************
112 : C VERSION : v1, 16/06/98, version initiale
113 : C HISTORIQUE : v1, 16/06/98, création
114 : C HISTORIQUE : 09/02/99 : ajout préconditionneur <> MATRIK
115 : C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau gmri0
116 : C (factorisations incomplètes stockées au format MSR (Modified Sparse
117 : C Row) (voir la doc de Sparskit version 2+)
118 : C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
119 : C HISTORIQUE : 09/04/04 : ajout des préconditionneurs ILUTP
120 : C HISTORIQUE :
121 : C HISTORIQUE :
122 : C***********************************************************************
123 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
124 : C en cas de modification de ce sous-programme afin de faciliter
125 : C la maintenance !
126 : C***********************************************************************
127 : *
128 : * .. Includes et pointeurs associés ..
129 : -INC CCOPTIO
130 : -INC SMLREEL
131 : INTEGER JG
132 : POINTEUR LRES.MLREEL
133 : -INC SMMATRIK
134 : POINTEUR MAPREC.MATRIK
135 : POINTEUR KMORS.PMORS
136 : POINTEUR KISA.IZA
137 : POINTEUR KS2B.IZA
138 : POINTEUR INCX.IZA
139 : POINTEUR INVDIA.IZA
140 : POINTEUR ILUM.PMORS
141 : POINTEUR ILUI.IZA
142 : * .. Parameters
143 : * This is a restart parameter for GMRES
144 : INTEGER RESTRT
145 : INTEGER LDW,NLDW
146 : SEGMENT SPACE
147 : REAL*8 WRK(LDW,NLDW)
148 : ENDSEGMENT
149 : * .. Work Vectors for GMRES..
150 : POINTEUR IV.SPACE
151 : POINTEUR IH.SPACE
152 : POINTEUR IR.IZA,IW.IZA,IAV.IZA
153 : POINTEUR IY.IZA,IGCOS.IZA,IGSIN.IZA
154 : POINTEUR IS.IZA
155 : * .. Scalar Arguments ..
156 : INTEGER ITER, KPREC, IMPR, IRET
157 : REAL*8 RESID
158 : INTEGER NBVA,NJA,NTT,NTTPRE
159 : * .. Executable Statements ..
160 : *
161 : IRET = 0
162 : *
163 : * On récupère les paramètres
164 : *
165 : SEGACT MATRIK
166 : SEGACT MAPREC
167 : IF (KSYM.EQ.0) THEN
168 : IF (IMPR.GT.0) THEN
169 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
170 : $ 'use a Conjugate Gradient instead !'
171 : ENDIF
172 : C IRET=-2
173 : C GOTO 9999
174 : ENDIF
175 : * Pour le préconditionneur
176 : ILUM =MAPREC.KIDMAT(6)
177 : ILUI =MAPREC.KIDMAT(7)
178 : IDMAT=MAPREC.KIDMAT(1)
179 : SEGACT IDMAT
180 : INVDIA=IDIAG
181 : SEGDES IDMAT
182 :
183 : SEGACT KMORS
184 : NTT =KMORS.IA(/1)-1
185 : NJA =KMORS.JA(/1)
186 : SEGACT KISA
187 : SEGACT KS2B
188 : SEGACT INCX*MOD
189 : C Paramètres des préconditionnements diagonaux et D-ILU
190 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
191 : C Est-il compatible avec la matrice ?
192 : SEGACT INVDIA
193 : NTTPRE=INVDIA.A(/1)
194 : IF (NTTPRE.NE.NTT) THEN
195 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
196 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
197 : WRITE(IOIMP,*) 'NTT=',NTT
198 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
199 : IRET=-2
200 : GOTO 9999
201 : ENDIF
202 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
203 : C ilutp, ilutpg, ilutpg2
204 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
205 : SEGACT ILUM
206 : SEGACT ILUI
207 : NTTPRE=ILUM.IA(/1)
208 : IF (NTTPRE.NE.NTT) THEN
209 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
210 : $ 'ne sont pas compatibles...'
211 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
212 : IRET=-2
213 : GOTO 9999
214 : ENDIF
215 : ELSE
216 : INVDIA=0
217 : ENDIF
218 : C
219 : C Initialisation de l'historique de convergence
220 : C
221 : JG=ITER
222 : SEGINI LRES
223 : C
224 : C
225 : C Initialisation des vecteurs de travail pour GMRES
226 : C
227 : NBVA=NTT
228 : SEGINI IR,IW,IAV
229 : NBVA=RESTRT
230 : SEGINI IY,IGCOS,IGSIN
231 : NBVA=RESTRT+1
232 : SEGINI IS
233 : LDW=NTT
234 : NLDW=RESTRT+1
235 : SEGINI IV
236 : LDW=RESTRT+1
237 : NLDW=RESTRT
238 : SEGINI IH
239 : C
240 : IF (KPREC.EQ.0) THEN
241 : C
242 : C Le GMRES non préconditionné (fortran pur)
243 : C
243 : CALL GMR(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
244 : $ KS2B.A,INCX.A,
245 : $ LRES.PROG,
246 : $ IR.A,IW.A,IAV.A,
247 : $ IY.A,IGCOS.A,IGSIN.A,IS.A,
248 : $ IV.WRK, IH.WRK,
249 : $ ITER,RESTRT,RESID,IMPR,IRET)
250 : ELSE IF (KPREC.EQ.1) THEN
251 : C
252 : C Le GMRES préconditionné par la diagonale (fortran pur)
253 : C
253 : CALL GMRD(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
254 : $ KS2B.A,INCX.A,
255 : $ INVDIA.A,
256 : $ LRES.PROG,
257 : $ IR.A,IW.A,IAV.A,
258 : $ IY.A,IGCOS.A,IGSIN.A,IS.A,
259 : $ IV.WRK, IH.WRK,
260 : $ ITER,RESTRT,RESID,IMPR,IRET)
261 : ELSE IF (KPREC.EQ.2) THEN
262 : C
263 : C Le GMRES préconditionné par une factorisation
264 : C incomplete D-ILU (fortran pur)
265 : C
265 : CALL GMRDI(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
266 : $ KS2B.A,INCX.A,
267 : $ INVDIA.A,
268 : $ LRES.PROG,
269 : $ IR.A,IW.A,IAV.A,
270 : $ IY.A,IGCOS.A,IGSIN.A,IS.A,
271 : $ IV.WRK, IH.WRK,
272 : $ ITER,RESTRT,RESID,IMPR,IRET)
273 : IF (IRET.NE.0) GOTO 9999
274 : ELSE IF (KPREC.GE.3.AND.KPREC.LE.9) THEN
275 : C
276 : C Le GMRES préconditionné par une factorisation
277 : C incomplete (fortran pur)
278 : C
278 : CALL GMRI0(NTT,NJA,KMORS.IA,KMORS.JA,KISA.A,
279 : $ KS2B.A,INCX.A,
280 : $ ILUM.IA,ILUM.JA,ILUI.A,
281 : $ LRES.PROG,
282 : $ IR.A,IW.A,IAV.A,
283 : $ IY.A,IGCOS.A,IGSIN.A,IS.A,
284 : $ IV.WRK, IH.WRK,
285 : $ ITER,RESTRT,RESID,IMPR,IRET)
286 : ELSE
287 : WRITE (IOIMP,*) 'KPREC out of range (0..9)'
288 : IRET=-4
289 : GOTO 9999
290 : ENDIF
291 : C
292 : C Désactivation-suppression
293 : C
294 : SEGSUP IH
295 : SEGSUP IV
296 : SEGSUP IS
297 : SEGSUP IY,IGCOS,IGSIN
298 : SEGSUP IR,IW,IAV
299 : JG=ITER
300 : SEGADJ LRES
301 : SEGDES LRES
302 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
303 : SEGDES INVDIA
304 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
305 : SEGDES ILUI
306 : SEGDES ILUM
307 : ENDIF
308 : SEGDES INCX
309 : SEGDES KS2B
310 : SEGDES KISA
311 : SEGDES KMORS
312 : SEGDES MAPREC
313 : SEGDES MATRIK
314 : C
315 : C A problem has occured in the GMRES method
316 : C
317 : IF (IRET.LT.0) GOTO 9999
318 : *
319 : * Normal termination
320 : *
321 : RETURN
322 : *
323 : * Format handling
324 : *
325 : 1002 FORMAT(10(1X,1PE11.4))
326 : *
327 : * Error handling
328 : *
329 : 9999 CONTINUE
330 : WRITE(IOIMP,*) 'An error was detected in prgmr.eso'
331 : RETURN
332 : *
333 : * End of prgmr
334 : *
335 : END
336 :
337 :
338 :
339 :
340 :
341 :
342 :
343 :
344 :
345 :
346 :
347 :
© Cast3M 2003 - All rights reserved.
Disclaimer