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