1 : C PRILTP SOURCE MAGN 09/07/30 21:20:31 6458
2 : SUBROUTINE PRILTP(AMORS,AISA,MATRIK,XLFIL,XDTOL,XSPIV,
3 : $ IVARI,
4 : $ INCX,NORMP,LRACOU,
5 : $ IMPR,IRET)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : IMPLICIT INTEGER (I-N)
8 : C***********************************************************************
9 : C NOM : PRILTP
10 : C DESCRIPTION :
11 : C Calcul du préconditionneur ILUTP d'une matrice Morse
12 : C Le préconditionneur est une matrice stockée
13 : C au format MSR (Modified Sparse Row, stockage de l'inverse de la
14 : C diagonale) de meme profil que la matrice Morse (format CSR) qu'il
15 : C préconditionne.
16 : C Le profil et les valeurs du préconditionneur sont
17 : C stockés dans KIDMAT(6 et 7) (réutilisation de l'existant).
18 : C
19 : C Ce sous-programme est en fait une interface à :
20 : C meiltp
21 : C qui est en Fortran presque pur (pour raison de rapidité)
22 : C et effectue la construction proprement dite du
23 : C préconditionneur.
24 : C
25 : C ATTENTION : pour une matrice A quelconque, la factorisation
26 : C --------- ILUTP peut ne pas exister (pivot nul) ou avoir
27 : C des pivots négatifs MEME SI la factorisation
28 : C complète de A existe et n'a que des pivots
29 : C positifs.
30 : C
31 : C
32 : C LANGAGE : ESOPE
33 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
34 : C mél : gounand@semt2.smts.cea.fr
35 : C REFERENCE (bibtex-like) :
36 : C @BOOK{templates,
37 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
38 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
39 : C H. Van der Vorst},
40 : C TITLE={Templates for the Solution of Linear Systems :
41 : C Building Blocks for Iterative Methods},
42 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
43 : C -> URL : http://www.netlib.org/templates/Templates.html
44 : C Sparskit : a basic tool kit for sparse matrix computations
45 : C Version 2 (Youcef Saad)
46 : C -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
47 : C REMARQUES :
48 : C On pourrait effectuer l'ordonnancement des colonnes du préconditionneur
49 : C (peut-être que la montée-descente va plus vite après, il faudrait
50 : C tester)
51 : C***********************************************************************
52 : C APPELES : MEILUT
53 : C APPELE PAR : KRES2
54 : C***********************************************************************
55 : C ENTREES : MATRIK, XLFIL, XDTOL, IMPR
56 : C SORTIES : ILUM, ILUI (KIDMAT(6-7) dans MATRIK), IRET
57 : C CODE RETOUR (IRET) : 0 si ok
58 : C <0 si problème
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 IMPR : niveau d'impression
63 : C***********************************************************************
64 : C VERSION : v1, 08/04/2004, version initiale
65 : C HISTORIQUE : v1, 08/04/2004, création
66 : C HISTORIQUE :
67 : C***********************************************************************
68 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
69 : C en cas de modification de ce sous-programme afin de faciliter
70 : C la maintenance !
71 : C***********************************************************************
72 : -INC CCOPTIO
73 : -INC CCREEL
74 : -INC SMMATRIK
75 : INTEGER NTT,NJA,NBVA
76 : POINTEUR AMORS.PMORS
77 : POINTEUR AISA.IZA
78 : POINTEUR KMORS.PMORS
79 : POINTEUR KISA.IZA
80 : POINTEUR ILUM.PMORS
81 : POINTEUR ILUI.IZA
82 : POINTEUR IDMATP.IDMAT
83 : POINTEUR IDMATD.IDMAT
84 : POINTEUR INCX.IZA
85 : POINTEUR INCX2.IZA
86 : -INC SMLENTI
87 : INTEGER JG
88 : POINTEUR IWORK.MLENTI
89 : POINTEUR JWORK.MLENTI
90 : POINTEUR KWORK.MLENTI
91 : POINTEUR IPERM.MLENTI
92 : POINTEUR JPERM.MLENTI
93 : -INC SMLREEL
94 : POINTEUR RWORK.MLREEL
95 : POINTEUR QWORK.MLREEL
96 : POINTEUR NORMP.MLREEL
97 : POINTEUR NORMP2.MLREEL
98 : * -INC SMLLOGI
99 : SEGMENT MLLOGI
100 : LOGICAL LOGI(JG)
101 : ENDSEGMENT
102 : POINTEUR LWORK.MLLOGI
103 : *
104 : REAL*8 XLFIL,XDTOL
105 : INTEGER IMPR,IRET
106 : *
107 : INTEGER ILFIL,NNZA,NNZLU,NNZMLU,NTTDDL
108 : REAL*8 XNTLIG
109 : LOGICAL LRACOU
110 : *
111 : IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans priltp.eso'
112 : C On récupère les segments utiles
113 : SEGACT MATRIK
114 : ILUM =KIDMAT(6)
115 : ILUI =KIDMAT(7)
116 : SEGDES MATRIK
117 : C Le préconditionneur est-il déjà construit ?
118 : IF ((ILUM.EQ.0).OR.(ILUI.EQ.0)) THEN
119 : C La matrice Morse et son préconditionneur ont le meme
120 : C profil
121 : *! Attention, on modifie la matrice dans meiltp
122 : SEGACT AMORS*MOD
123 : NTTDDL=AMORS.IA(/1)-1
124 : NNZA=AMORS.JA(/1)
125 : SEGACT AISA*MOD
126 : *
127 : * Ceci bugge avec ILUTPGO2 !!!!!
128 : *
129 : *
130 : * XNTLIG=DBLE(NNZA)/DBLE(NTTDDL)
131 : * ILFIL=INT(XNTLIG*XLFIL/2.D0)
132 : * NNZMLU=NTTDDL*(2*ILFIL+1)
133 : *
134 : NNZMLU = INT(DBLE(NNZA)*XLFIL) + 1
135 : *
136 : NTT=NTTDDL-1
137 : NJA=NNZMLU+1
138 : SEGINI ILUM
139 : NBVA=NNZMLU+1
140 : SEGINI ILUI
141 : JG=NTTDDL
142 : SEGINI IWORK
143 : JG=NTTDDL
144 : SEGINI JWORK
145 : JG=NTTDDL
146 : SEGINI RWORK
147 : JG=NTTDDL
148 : SEGINI IPERM
149 : JG=NTTDDL
150 : SEGINI JPERM
151 : * WRITE(IOIMP,*) 'NTTDDL=',NTTDDL,'NNZA=',NNZA,
152 : * $ 'NNZMLU=',NNZMLU
153 : * WRITE(IOIMP,*) 'ILFIL=',ILFIL,'XLFIL=',XLFIL,'XDTOL=',XDTOL
154 : * IF (XSPIV.LT.XZERO.OR.XSPIV.GT.1.D0) THEN
155 : * WRITE(IOIMP,*) 'Error : XSPIV=',XSPIV
156 : * GOTO 9999
157 : * ENDIF
158 : IF (IVARI.EQ.0) THEN
159 : C Les boucles sont en Fortran pur
159 : CALL MEILTP(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
160 : $ XLFIL,XDTOL,XSPIV,NNZMLU,
161 : $ ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
162 : $ IWORK.LECT,JWORK.LECT,RWORK.PROG,
163 : $ IPERM.LECT,JPERM.LECT,
164 : $ IMPR,IRET)
165 : IF (IRET.NE.0) GOTO 9999
166 : ELSEIF (IVARI.EQ.1) THEN
167 : JG=NTTDDL
168 : SEGINI LWORK
168 : CALL MILTP2(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
169 : $ XLFIL,XDTOL,XSPIV,NNZMLU,
170 : $ ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
171 : $ IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
172 : $ IPERM.LECT,JPERM.LECT,
173 : $ IMPR,IRET)
174 : IF (IRET.NE.0) GOTO 9999
175 : SEGSUP LWORK
176 : ELSEIF (IVARI.EQ.2) THEN
177 : JG=NTTDDL
178 : SEGINI LWORK
179 : JG=NTTDDL
180 : SEGINI KWORK
181 : JG=NTTDDL
182 : SEGINI QWORK
182 : CALL MILTP3(NTTDDL,NNZA,AISA.A,AMORS.JA,AMORS.IA,
183 : $ XLFIL,XDTOL,XSPIV,NNZMLU,
184 : $ ILUI.A,ILUM.JA,ILUM.IA,NNZLU,
185 : $ IWORK.LECT,JWORK.LECT,RWORK.PROG,LWORK.LOGI,
186 : $ KWORK.LECT,QWORK.PROG,
187 : $ IPERM.LECT,JPERM.LECT,
188 : $ IMPR,IRET)
189 : IF (IRET.NE.0) GOTO 9999
190 : SEGSUP QWORK
191 : SEGSUP KWORK
192 : SEGSUP LWORK
193 : ELSE
194 : WRITE(IOIMP,*) 'IVARI=',IVARI,' non prevu'
195 : ENDIF
196 : IF (IMPR.GT.1) THEN
197 : WRITE(IOIMP,*) 'Préconditionneur : nb.termesstockés=',
198 : $ NNZLU-1
199 : ENDIF
200 : *
201 : * SEGPRT,JPERM
202 : SEGSUP IPERM
203 : SEGSUP RWORK
204 : SEGSUP JWORK
205 : SEGSUP IWORK
206 : NBVA=NNZLU
207 : SEGADJ,ILUI
208 : SEGDES ILUI
209 : NTT=NTTDDL-1
210 : NJA=NNZLU
211 : SEGADJ,ILUM
212 : SEGDES ILUM
213 : SEGDES AISA
214 : SEGDES AMORS
215 : * SEGPRT,ILUM
216 : * SEGPRT,ILUI
217 : *
218 : * On s'occupe des permutations
219 : *
220 : SEGACT MATRIK*MOD
221 : IF (.NOT.LRACOU) THEN
222 : KMORS=MATRIK.KIDMAT(4)
223 : SEGACT KMORS*MOD
224 : DO IJA=KMORS.IA(1),KMORS.IA(NTTDDL+1)-1
225 : KMORS.JA(IJA)=JPERM.LECT(KMORS.JA(IJA))
226 : ENDDO
227 : SEGDES KMORS
228 : ENDIF
229 : IDMATP=MATRIK.KIDMAT(1)
230 : IDMATD=MATRIK.KIDMAT(2)
231 : IF (IDMATD.EQ.0) THEN
232 : WRITE(IOIMP,*) 'Erreur hyper grave'
233 : ENDIF
234 : IF (IDMATP.EQ.IDMATD) THEN
235 : SEGINI,IDMATP=IDMATD
236 : ELSE
237 : SEGACT IDMATP*MOD
238 : ENDIF
239 : * Copie de NUAN dans NUNA
240 : DO ITTDDL=1,NTTDDL
241 : IDMATP.NUNA(ITTDDL)=IDMATP.NUAN(ITTDDL)
242 : ENDDO
243 : * Mise à jour de NUAN
244 : DO ITTDDL=1,NTTDDL
245 : IDMATP.NUAN(ITTDDL)=JPERM.LECT(IDMATP.NUNA(ITTDDL))
246 : ENDDO
247 : * Mise à jour de NUNA
248 : DO ITTDDL=1,NTTDDL
249 : IDMATP.NUNA(IDMATP.NUAN(ITTDDL))=ITTDDL
250 : ENDDO
251 : INCX2=INCX
252 : SEGACT INCX2
253 : NBVA=NTTDDL
254 : SEGINI,INCX
255 : DO ITTDDL=1,NTTDDL
256 : INCX.A(JPERM.LECT(ITTDDL))=INCX2.A(ITTDDL)
257 : ENDDO
258 : SEGDES INCX
259 : SEGSUP INCX2
260 : IF (NORMP.NE.0) THEN
261 : NORMP2=NORMP
262 : SEGACT NORMP2
263 : JG=NTTDDL
264 : SEGINI,NORMP
265 : DO ITTDDL=1,NTTDDL
266 : NORMP.PROG(JPERM.LECT(ITTDDL))=NORMP2.PROG(ITTDDL)
267 : ENDDO
268 : SEGDES NORMP
269 : SEGSUP NORMP2
270 : ENDIF
271 : SEGSUP JPERM
272 : C SEGDES IDMATP
273 : SEGDES IDMATP
274 : MATRIK.KIDMAT(1)=IDMATP
275 : SEGDES MATRIK
276 : C
277 : C On stocke la factorisation obtenue du préconditionneur
278 : C
279 : SEGACT MATRIK*MOD
280 : KIDMAT(6)=ILUM
281 : KIDMAT(7)=ILUI
282 : SEGDES MATRIK
283 : C
284 : IF (IMPR.GT.6) THEN
285 : WRITE(IOIMP,*) 'création du préconditionneur Morse',
286 : $ ' de pointeurs',ILUM,'et',ILUI
287 : IF (IMPR.GT.8) THEN
287 : CALL ECMORS(ILUM,ILUI,(IMPR-1))
288 : ENDIF
289 : ENDIF
290 : ELSE
291 : IF (IMPR.GT.1) THEN
292 : WRITE(IOIMP,*)
293 : $ 'Préconditionneur deja construit.'
294 : ENDIF
295 : IF (IMPR.GT.6) THEN
296 : WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :',
297 : $ 'ILUM=',ILUM,' et ILUI=',ILUI
298 : IF (IMPR.GT.8) THEN
298 : CALL ECMORS(ILUM,ILUI,(IMPR-1))
299 : ENDIF
300 : ENDIF
301 : ENDIF
302 : *
303 : * Normal termination
304 : *
305 : RETURN
306 : *
307 : * Format handling
308 : *
309 : *
310 : * Error handling
311 : *
312 : 9999 CONTINUE
313 : WRITE(IOIMP,*) 'An error was detected in priltp.eso'
314 : IRET=1
315 : RETURN
316 : *
317 : * End of PRILTP
318 : *
319 : END
320 :
321 :
322 :
323 :
324 :
325 :
326 :
327 :
328 :
© Cast3M 2003 - All rights reserved.
Disclaimer