1 : C PRGMR2 SOURCE GOUNAND 11/10/07 21:16:09 7149
2 : SUBROUTINE PRGMR2(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,INCX,
3 : $ ITER,INMV,RESID,KPREC,
4 : $ RESTRT,ICALRS,IDDOT,IMVEC,IMPR,IRET)
5 : IMPLICIT INTEGER(I-N)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : C***********************************************************************
8 : C NOM : PRGMR2
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 : 10/02/06 : Nettoyage
121 : C HISTORIQUE :
122 : C HISTORIQUE :
123 : C***********************************************************************
124 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
125 : C en cas de modification de ce sous-programme afin de faciliter
126 : C la maintenance !
127 : C***********************************************************************
128 : *
129 : * .. Includes et pointeurs associés ..
130 : -INC CCOPTIO
131 : -INC SMLREEL
132 : INTEGER JG
133 : POINTEUR LRES.MLREEL
134 : -INC SMLENTI
135 : POINTEUR LNMV.MLENTI
136 : -INC SMMATRIK
137 : POINTEUR MAPREC.MATRIK
138 : POINTEUR KMORS.PMORS
139 : POINTEUR KISA.IZA
140 : POINTEUR KS2B.IZA
141 : POINTEUR INCX.IZA
142 : POINTEUR INVDIA.IZA
143 : POINTEUR ILUM.PMORS
144 : POINTEUR ILUI.IZA
145 : * .. Parameters
146 : * This is a restart parameter for GMRES
147 : INTEGER RESTRT
148 : INTEGER LDW,NLDW
149 : SEGMENT SPACE
150 : REAL*8 WRK(LDW,NLDW)
151 : ENDSEGMENT
152 : SEGMENT SPACE2
153 : POINTEUR WRK2(NIZA).IZA
154 : ENDSEGMENT
155 : * .. Work Vectors for GMRES..
156 : POINTEUR ITMP.IZA
157 : POINTEUR IV.SPACE2
158 : POINTEUR IH.SPACE
159 : POINTEUR IR.IZA,IW.IZA,IAV.IZA
160 : POINTEUR IY.IZA,IGCOS.IZA,IGSIN.IZA
161 : POINTEUR IS.IZA
162 : * .. Scalar Arguments ..
163 : INTEGER ITER, KPREC, IMPR, IRET
164 : REAL*8 RESID
165 : INTEGER NBVA,NJA,NTT,NTTPRE
166 : * .. Executable Statements ..
167 : *
168 : IRET = 0
169 : *
170 : * On récupère les paramètres
171 : *
172 : SEGACT MATRIK
173 : SEGACT MAPREC
174 : IF (KSYM.EQ.0) THEN
175 : IF (IMPR.GT.0) THEN
176 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
177 : $ 'use a Conjugate Gradient instead !'
178 : ENDIF
179 : C IRET=-2
180 : C GOTO 9999
181 : ENDIF
182 : * Pour le préconditionneur
183 : ILUM =MAPREC.KIDMAT(6)
184 : ILUI =MAPREC.KIDMAT(7)
185 : IDMAT=MAPREC.KIDMAT(1)
186 : SEGACT IDMAT
187 : INVDIA=IDIAG
188 : SEGDES IDMAT
189 :
190 : SEGACT KMORS
191 : NTT =KMORS.IA(/1)-1
192 : * NJA =KMORS.JA(/1)
193 : SEGACT KISA
194 : SEGACT KS2B
195 : SEGACT INCX*MOD
196 : C Paramètres des préconditionnements diagonaux et D-ILU
197 : IF (KPREC.NE.0) THEN
198 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
199 : C Est-il compatible avec la matrice ?
200 : SEGACT INVDIA
201 : NTTPRE=INVDIA.A(/1)
202 : IF (NTTPRE.NE.NTT) THEN
203 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
204 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
205 : WRITE(IOIMP,*) 'NTT=',NTT
206 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
207 : IRET=-2
208 : GOTO 9999
209 : ENDIF
210 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
211 : C ilutp, ilutpg, ilutpg2
212 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
213 : SEGACT ILUM
214 : SEGACT ILUI
215 : NTTPRE=ILUM.IA(/1)
216 : IF (NTTPRE.NE.NTT) THEN
217 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
218 : $ 'ne sont pas compatibles...'
219 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
220 : IRET=-2
221 : GOTO 9999
222 : ENDIF
223 : ENDIF
224 : ENDIF
225 : C
226 : C Initialisation de l'historique de convergence
227 : C
228 : JG=ITER+1
229 : SEGINI LNMV
230 : SEGINI LRES
231 : C
232 : C
233 : C Initialisation des vecteurs de travail pour GMRES
234 : C
235 : NBVA=NTT
236 : SEGINI IR,IW,IAV
237 : NBVA=RESTRT
238 : SEGINI IY,IGCOS,IGSIN
239 : NBVA=RESTRT+1
240 : SEGINI IS
241 : NBVA=NTT
242 : NIZA=RESTRT+1
243 : SEGINI IV
244 : DO I=1,NIZA
245 : SEGINI ITMP
246 : IV.WRK2(I)=ITMP
247 : ENDDO
248 : LDW=RESTRT+1
249 : NLDW=RESTRT
250 : SEGINI IH
251 : C
251 : CALL GMR2(KMORS,KISA,KS2B,INCX,
252 : $ KPREC,INVDIA,ILUM,ILUI,
253 : $ LRES,LNMV,
254 : $ IR,IW,IAV,
255 : $ IY,IGCOS,IGSIN,IS,
256 : $ IV, IH,
257 : $ ITER,INMV,RESTRT,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
258 : C
259 : C Désactivation-suppression
260 : C
261 : SEGSUP IH
262 : DO I=1,NIZA
263 : ITMP=IV.WRK2(I)
264 : SEGSUP ITMP
265 : ENDDO
266 : SEGSUP IV
267 : SEGSUP IS
268 : SEGSUP IY,IGCOS,IGSIN
269 : SEGSUP IR,IW,IAV
270 : JG=ITER+1
271 : SEGADJ LRES
272 : SEGDES LRES
273 : SEGADJ LNMV
274 : SEGDES LNMV
275 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
276 : SEGDES INVDIA
277 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
278 : SEGDES ILUI
279 : SEGDES ILUM
280 : ENDIF
281 : SEGDES INCX
282 : SEGDES KS2B
283 : SEGDES KISA
284 : SEGDES KMORS
285 : SEGDES MAPREC
286 : SEGDES MATRIK
287 : C
288 : C A problem has occured in the GMRES method
289 : C
290 : IF (IRET.LT.0) GOTO 9999
291 : *
292 : * Normal termination
293 : *
294 : RETURN
295 : *
296 : * Format handling
297 : *
298 : *
299 : * Error handling
300 : *
301 : 9999 CONTINUE
302 : WRITE(IOIMP,*) 'An error was detected in prgmr2.eso'
303 : RETURN
304 : *
305 : * End of prgmr2
306 : *
307 : END
308 :
309 :
310 :
311 :
312 :
313 :
314 :
315 :
316 :
317 :
318 :
319 :
320 :
321 :
© Cast3M 2003 - All rights reserved.
Disclaimer