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