1 : C PRBCG SOURCE GOUNAND 11/10/07 21:16:04 7149
2 : SUBROUTINE PRBCG(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 : PRBCG
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 : POINTEUR IRTLD.IZA,IPTLD.IZA,IQTLD.IZA,IZTLD.IZA
131 : * .. Scalar Arguments ..
132 : INTEGER ITER, KPREC, IMPR, IRET
133 : REAL*8 RESID
134 : INTEGER NBVA,NJA,NTT,NTTPRE
135 : * .. Executable Statements ..
136 : *
137 : IRET = 0
138 : *
139 : * On récupère les paramètres
140 : *
141 : SEGACT MATRIK
142 : SEGACT MAPREC
143 : IF (KSYM.NE.0) THEN
144 : IF (IMPR.GT.0) THEN
145 : WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
146 : $ 'use Bi-CGSTAB instead !'
147 : * ENDIF
148 : * IRET=-1
149 : * GOTO 9999
150 : ENDIF
151 : ENDIF
152 : * Pour le préconditionneur
153 : ILUM =MAPREC.KIDMAT(6)
154 : ILUI =MAPREC.KIDMAT(7)
155 : IDMAT=MAPREC.KIDMAT(1)
156 : SEGACT IDMAT
157 : INVDIA=IDIAG
158 : SEGDES IDMAT
159 : SEGACT KMORS
160 : NTT =KMORS.IA(/1)-1
161 : * NJA =KMORS.JA(/1)
162 : SEGACT KISA
163 : SEGACT KS2B
164 : SEGACT INCX*MOD
165 : C Paramètres des préconditionnements diagonaux et D-ILU
166 : IF (KPREC.NE.0) THEN
167 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
168 : C Est-il compatible avec la matrice ?
169 : SEGACT INVDIA
170 : NTTPRE=INVDIA.A(/1)
171 : IF (NTTPRE.NE.NTT) THEN
172 : WRITE(IOIMP,*) 'La matrice et son préconditionnement'
173 : WRITE(IOIMP,*) 'ne sont pas compatibles...'
174 : WRITE(IOIMP,*) 'NTT=',NTT
175 : WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
176 : IRET=-2
177 : GOTO 9999
178 : ENDIF
179 : C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
180 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
181 : SEGACT ILUM
182 : SEGACT ILUI
183 : NTTPRE=ILUM.IA(/1)
184 : IF (NTTPRE.NE.NTT) THEN
185 : WRITE(IOIMP,*) 'La matrice et son préconditionnement',
186 : $ 'ne sont pas compatibles...'
187 : WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
188 : IRET=-2
189 : GOTO 9999
190 : ENDIF
191 : ELSE
192 : WRITE(IOIMP,*) 'Erreur de programmation'
193 : GOTO 9999
194 : ENDIF
195 : ENDIF
196 : C
197 : C Initialisation de l'historique de convergence
198 : C
199 : JG=ITER+1
200 : SEGINI LNMV
201 : SEGINI LRES
202 : C
203 : C
204 : C Initialisation des vecteurs de travail pour le gradient conjugué
205 : C
206 : NBVA=NTT
207 : SEGINI IR,IRTLD,IZ,IZTLD,IP,IPTLD
208 : IQ=IZ
209 : IQTLD=IZTLD
210 : C
211 : C Appel du gradient conjugué
212 : C
212 : CALL BICG(KMORS,KISA,KS2B,INCX,
213 : $ KPREC,INVDIA,ILUM,ILUI,
214 : $ LRES,LNMV,
215 : $ IR,IRTLD,IZ,IZTLD,IP,IPTLD,IQ,IQTLD,
216 : $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
217 : C
218 : C Désactivation
219 : C
220 : SEGSUP IR,IRTLD,IZ,IZTLD,IP,IPTLD
221 : JG=ITER+1
222 : SEGADJ LRES
223 : SEGDES LRES
224 : SEGADJ LNMV
225 : SEGDES LNMV
226 : IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
227 : SEGDES INVDIA
228 : ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
229 : SEGDES ILUI
230 : SEGDES ILUM
231 : ENDIF
232 : SEGDES INCX
233 : SEGDES KS2B
234 : SEGDES KISA
235 : SEGDES KMORS
236 : SEGDES MAPREC
237 : SEGDES MATRIK
238 : *
239 : * Normal termination
240 : *
241 : RETURN
242 : *
243 : * Format handling
244 : *
245 : *
246 : * Error handling
247 : *
248 : 9999 CONTINUE
249 : WRITE(IOIMP,*) 'An error was detected in prbcg.eso'
250 : RETURN
251 : *
252 : * End of PRBCG.
253 : *
254 : END
255 :
256 :
257 :
258 :
259 :
260 :
261 :
262 :
263 :
264 :
265 :
266 :
267 :
268 :
269 :
270 :
© Cast3M 2003 - All rights reserved.
Disclaimer