1 : C MEILUT SOURCE GOUNAND 05/06/06 21:15:39 5078
2 : SUBROUTINE MEILUT(NTT,NNZA,A,JA,IA,
3 : $ XLFIL,XDTOL,NNZMLU,
4 : $ ALU,JLU,JU,NNZLU,
5 : $ IWORK,JWORK,RWORK,
6 : $ IVARI,IMPR,IRET)
7 : IMPLICIT REAL*8 (A-H,O-Z)
8 : IMPLICIT INTEGER (I-N)
9 : C***********************************************************************
10 : C NOM : MEILUT
11 : C DESCRIPTION :
12 : C Calcul du préconditionneur ILUT d'une matrice Morse (IVARI=0)
13 : C ou d'une variante qui remplit mieux la mémoire et fonctionne mieux
14 : C quelquefois (IVARI=1).
15 : C ILUT : Incomplete LU factorization with dual truncation strategy
16 : C La matrice Morse est au format CSR (Column Sparse Row)
17 : C Le préconditionneur est au format MSR (Modified Sparse Row,
18 : C stockage de l'inverse de la diagonale) (ALU, JLU).
19 : C JU est un tableau supplémentaire pour le repérage des diagonales.
20 : C (voir la doc. de Sparskit version 2)
21 : C
22 : C LANGAGE : FORTRAN 77 (sauf E/S)
23 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
24 : C mél : gounand@semt2.smts.cea.fr
25 : C REFERENCE (bibtex-like) :
26 : C @BOOK{templates,
27 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
28 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
29 : C H. Van der Vorst},
30 : C TITLE={Templates for the Solution of Linear Systems :
31 : C Building Blocks for Iterative Methods},
32 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
33 : C -> URL : http://www.netlib.org/templates/Templates.html
34 : C Sparskit : a basic tool kit for sparse matrix computations
35 : C Version 2 (Youcef Saad)
36 : C -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
37 : C REMARQUES :
38 : C Quelques parties pourraient être améliorés :
39 : C - l'algorithme de quicksplit (cf.qsplit.eso)
40 : C - la partie où l'on doit effectuer l'élimination dans le bon ordre
41 : C pourrait utiliser un algorithme de tri type quicksort
42 : C
43 : C***********************************************************************
44 : C APPELES : QSPLIT
45 : C APPELE PAR : PRILUT
46 : C***********************************************************************
47 : C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
48 : C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
49 : C SORTIES : ALU,JLU,JU,NNZLU
50 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
51 : C***********************************************************************
52 : C VERSION : v1, 20/03/00, réécriture totale
53 : C VERSION : v1, 23/02/00, version initiale
54 : C HISTORIQUE : v1, 23/02/00, création
55 : C HISTORIQUE :
56 : C HISTORIQUE :
57 : C***********************************************************************
58 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
59 : C en cas de modification de ce sous-programme afin de faciliter
60 : C la maintenance !
61 : C***********************************************************************
62 : -INC CCOPTIO
63 : -INC CCREEL
64 : * .. Scalar Arguments ..
65 : INTEGER NTT,NNZA,NNZMLU,NNZLU
66 : INTEGER IVARI,IMPR,IRET
67 : REAL*8 XDTOL,XLFIL
68 : * ..
69 : * .. Array Arguments ..
70 : REAL*8 A(NNZA),ALU(NNZMLU+1)
71 : INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
72 : INTEGER IWORK(NTT)
73 : INTEGER JWORK(NTT)
74 : REAL*8 RWORK(NTT)
75 : *
76 : * .. Variables locales
77 : * .. Parameters
78 : REAL*8 ZERO ,ONE
79 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
80 : * ..
81 : C Nombre de pivots petits
82 : INTEGER NBPIVP
83 : C Nombre de pivots inférieurs à 0
84 : INTEGER NBPIVI
85 : *
86 : INTEGER ILFIL,ILFILL,ILFILU
87 : INTEGER COL,COLLU
88 : INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
89 : INTEGER JJL,JPOS,JROW,KJL,KLU
90 : INTEGER LEN,LENL,LENU
91 : INTEGER JU0,ITMP
92 : REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
93 : REAL*8 XNTLIG
94 : * REAL*8 XPOND1,XPOND2
95 : *
96 : * Executable statements
97 : *
98 : IF (XLFIL.LT.0.D0) THEN
99 : WRITE(IOIMP,*) 'XLFIL.LT.0.D0 non acceptable'
100 : GOTO 9999
101 : ENDIF
102 : PETIT=XZPREC**(0.75D0)
103 : NBPIVP=0
104 : NBPIVI=0
105 : c-----------------------------------------------------------------------
106 : c initialize ju0 (points to next element to be added to alu,jlu)
107 : c and pointer array.
108 : c-----------------------------------------------------------------------
109 : JU0=NTT+2
110 : JLU(1)=JU0
111 : *
112 : * On boucle sur les ddl
113 : *
114 : DO 1 ITT=1,NTT
115 : * Calcul de la norme de la ittème ligne
116 : ISTRT=IA(ITT)
117 : ISTOP=IA(ITT+1)-1
118 : TNORM=ZERO
119 : DO 11 IJA=ISTRT,ISTOP
120 : TNORM=TNORM+ABS(A(IJA))
121 : 11 CONTINUE
122 : TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
123 : IF (TNORM.LE.XPETIT) THEN
124 : WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
125 : GOTO 9999
126 : ENDIF
127 : * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
128 : * parties triangulaires inférieures et supérieurs.
129 : *
130 : * On peut le programmer autrement compte tenu du fait que les
131 : * colonnes sont ordonnées
132 : *
133 : LENL=0
134 : LENU=0
135 : DO 12 IJA=ISTRT,ISTOP
136 : COL=JA(IJA)
137 : VAL=A(IJA)
138 : IF (COL.LT.ITT) THEN
139 : LENL=LENL+1
140 : IWORK(LENL)=COL
141 : RWORK(LENL)=VAL
142 : JWORK(COL) =LENL
143 : ELSEIF (COL.EQ.ITT) THEN
144 : IWORK(ITT)=ITT
145 : RWORK(ITT)=VAL
146 : JWORK(ITT)=ITT
147 : ELSE
148 : LENU=LENU+1
149 : IWORK(ITT+LENU)=COL
150 : RWORK(ITT+LENU)=VAL
151 : JWORK(COL) =ITT+LENU
152 : ENDIF
153 : 12 CONTINUE
154 : IF (IVARI.EQ.1) THEN
155 : * Stratégie : on remplit autant que possible
156 : * normale meilug
157 : XNTLIG=DBLE(NNZMLU-JU0+1)/DBLE(NTT-ITT+1)
158 : * on retire la diagonale
159 : * ILFIL=INT(XNTLIG/2.D0)
160 : ILFIL=INT((XNTLIG-1.D0)/2.D0)
161 : * Stratégie meilug pondérée globalement
162 : * XNTLIG=(DBLE(NNZMLU-JU0+1)/DBLE(NTT-ITT+1))
163 : * $ *(1.5D0-0.5D0*(DBLE(ITT-1)/DBLE(NTT-1)))
164 : * ILFIL=INT(XNTLIG/2.D0)
165 : * Stratégie normale (idem meilut)
166 : ELSEIF (IVARI.EQ.0) THEN
167 : XNTLIG=DBLE(NNZA)/DBLE(NTT)
168 : * Pour faire comme dans meiltp
169 : * ILFIL=INT(XNTLIG*XLFIL/2.D0)
170 : ILFIL=INT((XNTLIG-1.D0)*XLFIL/2.D0)
171 : ELSE
172 : WRITE(IOIMP,*) 'Erreur de programmation IVARI=',IVARI
173 : ENDIF
174 : * On a enlevé les :
175 : * Corrections pour les extrémités de la matrice
176 : * IF (ILFIL.GT.ITT-1) THEN
177 : * ILFILL=ITT-1
178 : * ILFILU=(2*ILFIL)-ILFILL
179 : * ELSEIF (ILFIL.GT.(NTT-ITT)) THEN
180 : * ILFILU=NTT-ITT
181 : * ILFILL=(2*ILFIL)-ILFILU
182 : * ELSE
183 : * ILFILL=ILFIL
184 : * ILFILU=ILFIL
185 : * ENDIF
186 : * WRITE(IOIMP,*) 'itt,ilfil,ilfill,ilfilu=',ITT,ILFIL,ILFILL
187 : * $ ,ILFILU
188 : * Stratégie normale
189 : ILFILL=ILFIL
190 : ILFILU=ILFIL
191 : * Stratégie pondérée localement
192 : * XPOND1=DBLE(LENL)/DBLE(LENL+LENU)
193 : * XPOND2=DBLE(LENU)/DBLE(LENL+LENU)
194 : * ILFILL=INT(XPOND1*DBLE(ILFIL*2))
195 : * ILFILU=(ILFIL*2)-ILFILL
196 : * Elimination avec les lignes précédentes
197 : JJL=0
198 : LEN=0
199 : 13 CONTINUE
200 : JJL=JJL+1
201 : IF (JJL.LE.LENL) THEN
202 : c-----------------------------------------------------------------------
203 : c in order to do the elimination in the correct order we must select
204 : c the smallest column index among jw(k), k=jj+1, ..., lenl.
205 : c-----------------------------------------------------------------------
206 : JROW=IWORK(JJL)
207 : KJL=JJL
208 : c determine smallest column index
209 : DO 132 ILENL=JJL+1,LENL
210 : IF (IWORK(ILENL).LT.JROW) THEN
211 : JROW=IWORK(ILENL)
212 : KJL=ILENL
213 : ENDIF
214 : 132 CONTINUE
215 : IF (KJL.NE.JJL) THEN
216 : ITMP=IWORK(JJL)
217 : IWORK(JJL)=IWORK(KJL)
218 : IWORK(KJL)=ITMP
219 : JWORK(JROW)=JJL
220 : JWORK(ITMP)=KJL
221 : RTMP=RWORK(JJL)
222 : RWORK(JJL)=RWORK(KJL)
223 : RWORK(KJL)=RTMP
224 : ENDIF
225 : c zero out element in row by setting jw(n+jrow) to zero.
226 : JWORK(JROW)=0
227 : c get the multiplier for row to be eliminated (jrow).
228 : FACT=RWORK(JJL)*ALU(JROW)
229 : IF (ABS(FACT).LE.XDTOL) THEN
230 : GOTO 13
231 : ENDIF
232 : c
233 : c combine current row and row jrow
234 : c
235 : DO 134 KLU=JU(JROW),JLU(JROW+1)-1
236 : TERMLU=FACT*ALU(KLU)
237 : COLLU =JLU(KLU)
238 : JPOS =JWORK(COLLU)
239 : IF (COLLU.GE.ITT) THEN
240 : c dealing with upper part (including diagonal)
241 : IF (JPOS.EQ.0) THEN
242 : c this is a fill-in element
243 : LENU=LENU+1
244 : IF (LENU.GT.NTT-ITT) THEN
245 : WRITE(IOIMP,*) 'pb. matrice U'
246 : GOTO 9999
247 : ENDIF
248 : IWORK(ITT+LENU)=COLLU
249 : JWORK(COLLU)=ITT+LENU
250 : RWORK(ITT+LENU)=-TERMLU
251 : ELSE
252 : c this is not a fill-in element
253 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
254 : ENDIF
255 : ELSE
256 : c dealing with lower part.
257 : IF (JPOS.EQ.0) THEN
258 : c this is a fill-in element
259 : LENL=LENL+1
260 : IF (LENL.GT.ITT-1) THEN
261 : WRITE(IOIMP,*) 'pb. matrice L'
262 : GOTO 9999
263 : ENDIF
264 : IWORK(LENL) =COLLU
265 : JWORK(COLLU)=LENL
266 : RWORK(LENL) =-TERMLU
267 : ELSE
268 : c this is not a fill-in element
269 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
270 : ENDIF
271 : ENDIF
272 : 134 CONTINUE
273 : c store this pivot element -- (from left to right -- no danger of
274 : c overlap with the working elements in L (pivots).
275 : LEN=LEN+1
276 : IWORK(LEN)=JROW
277 : RWORK(LEN)=FACT
278 : GOTO 13
279 : ENDIF
280 : c reset double-pointer to zero (U-part including diagonal)
281 : DO 14 ILENU=0,LENU
282 : JWORK(IWORK(ITT+ILENU))=0
283 : 14 CONTINUE
284 : c update L-matrix
285 : LENL=LEN
286 : LEN=MIN(LENL,ILFILL)
287 : c sort by quick-split
287 : CALL QSPLIT(RWORK,IWORK,LENL,LEN,
288 : $ IMPR,IRET)
289 : c store L-part
290 : DO 15 ILEN=1,LEN
291 : JLU(JU0)=IWORK(ILEN)
292 : ALU(JU0)=RWORK(ILEN)
293 : JU0=JU0+1
294 : 15 CONTINUE
295 : c save pointer to beginning of row ii of U
296 : JU(ITT)=JU0
297 : c update U-matrix -- first apply dropping strategy
298 : LEN=0
299 : DO 16 ILEN=1,LENU
300 : IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
301 : LEN=LEN+1
302 : IWORK(ITT+LEN)=IWORK(ITT+ILEN)
303 : RWORK(ITT+LEN)=RWORK(ITT+ILEN)
304 : ENDIF
305 : 16 CONTINUE
306 : LENU=LEN
307 : LEN=MIN(LENU,ILFILU)
308 : c sort by quick-split
308 : CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU,LEN,
309 : $ IMPR,IRET)
310 : c copy
311 : DO 17 ILEN=ITT+1,ITT+LEN
312 : JLU(JU0)=IWORK(ILEN)
313 : ALU(JU0)=RWORK(ILEN)
314 : JU0=JU0+1
315 : 17 CONTINUE
316 : c store inverse of diagonal element of u
317 : VALPIV=RWORK(ITT)
318 : IF (VALPIV.EQ.ZERO) THEN
319 : WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
320 : $ 'le préconditionnement ILUT est impossible.'
321 : IRET=-2
322 : GOTO 9999
323 : ELSE
324 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
325 : NBPIVP=NBPIVP+1
326 : ENDIF
327 : IF (VALPIV.LT.ZERO) THEN
328 : NBPIVI=NBPIVI+1
329 : ENDIF
330 : ENDIF
331 : ALU(ITT)=ONE/RWORK(ITT)
332 : c update pointer to beginning of next row of U.
333 : JLU(ITT+1)=JU0
334 : 1 CONTINUE
335 : NNZLU=JLU(NTT+1)-1
336 : *
337 : * Warning(s)
338 : *
339 : IF (IMPR.GT.0) THEN
340 : IF (NBPIVP.GT.0) THEN
341 : WRITE(IOIMP,*) 'WARNING !',
342 : $ ' meilut :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
343 : ENDIF
344 : *! IF (NBPIVI.GT.0) THEN
345 : *! WRITE(IOIMP,*) 'WARNING !',
346 : *! $ ' meilut :',NBPIVI,' pivot(s) < 0'
347 : *! ENDIF
348 : ENDIF
349 : *
350 : * Normal termination
351 : *
352 : IRET=0
353 : RETURN
354 : *
355 : * Error handling
356 : *
357 : 9999 CONTINUE
358 : IRET=1
359 : WRITE(IOIMP,*) 'An error was detected in meilut.eso'
360 : RETURN
361 : *
362 : * End of MEILUT
363 : *
364 : END
365 :
366 :
367 :
368 :
369 :
370 :
371 :
© Cast3M 2003 - All rights reserved.
Disclaimer