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