1 : C MILTP3 SOURCE GOUNAND 05/02/16 21:16:19 5029
2 : SUBROUTINE MILTP3(NTT,NNZA,A,JA,IA,
3 : $ XLFIL,XDTOL,XSPIV,NNZMLU,
4 : $ ALU,JLU,JU,NNZLU,
5 : $ IWORK,JWORK,RWORK,LWORK,
6 : $ KWORK,QWORK,
7 : $ IPERM,JPERM,
8 : $ IMPR,IRET)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : IMPLICIT INTEGER (I-N)
11 : C***********************************************************************
12 : C NOM : MILTP3
13 : C DESCRIPTION :
14 : C Calcul du préconditionneur ILUTP d'une matrice Morse
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 On garde absolument tous les termes qui seraient dans ILU(0).
45 : C
46 : C
47 : C***********************************************************************
48 : C APPELES : QSPLIT
49 : C APPELE PAR : PRILUT
50 : C***********************************************************************
51 : C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
52 : C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
53 : C SORTIES : ALU,JLU,JU,NNZLU
54 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
55 : C***********************************************************************
56 : C VERSION : v1, 26/04/04, version initiale
57 : C HISTORIQUE : v1, 26/04/04, création
58 : C HISTORIQUE :
59 : C HISTORIQUE :
60 : C***********************************************************************
61 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
62 : C en cas de modification de ce sous-programme afin de faciliter
63 : C la maintenance !
64 : C***********************************************************************
65 : -INC CCOPTIO
66 : -INC CCREEL
67 : * .. Scalar Arguments ..
68 : INTEGER NTT,NNZA,NNZMLU,NNZLU
69 : INTEGER IMPR,IRET
70 : REAL*8 XDTOL,XLFIL
71 : * ..
72 : * .. Array Arguments ..
73 : REAL*8 A(NNZA),ALU(NNZMLU+1)
74 : INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
75 : INTEGER IWORK(NTT)
76 : INTEGER JWORK(NTT)
77 : REAL*8 RWORK(NTT)
78 : INTEGER KWORK(NTT)
79 : REAL*8 QWORK(NTT)
80 : LOGICAL LWORK(NTT)
81 : INTEGER IPERM(NTT)
82 : INTEGER JPERM(NTT)
83 : *
84 : * .. Variables locales
85 : * .. Parameters
86 : REAL*8 ZERO ,ONE
87 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
88 : * ..
89 : C Nombre de pivots petits
90 : INTEGER NBPIVP
91 : C Nombre de pivots inférieurs à 0
92 : INTEGER NBPIVI
93 : *
94 : INTEGER ILFIL,ILFILL,ILFILU
95 : INTEGER COL,COLLU
96 : INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
97 : INTEGER JJL,JPOS,JROW,KJL,KLU
98 : INTEGER LEN,LENL,LENU
99 : INTEGER JU0,ITMP
100 : REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
101 : LOGICAL LTMP,LFACT,LPIV,LFOUND
102 : REAL*8 XNTLIG
103 : * REAL*8 XPOND1,XPOND2
104 : *
105 : * Executable statements
106 : *
107 : IF (XLFIL.LT.1.D0) THEN
108 : WRITE(IOIMP,*) 'XLFIL.LT.1.D0 non acceptable'
109 : GOTO 9999
110 : ENDIF
111 : PETIT=XZPREC**(0.75D0)
112 : NBPIVP=0
113 : NBPIVI=0
114 : c-----------------------------------------------------------------------
115 : c initialize ju0 (points to next element to be added to alu,jlu)
116 : c and pointer array.
117 : c-----------------------------------------------------------------------
118 : JU0=NTT+2
119 : JLU(1)=JU0
120 : *
121 : * Initialize permutation arrays and lwork
122 : *
123 : DO ITT=1,NTT
124 : JWORK(ITT)=0
125 : IPERM(ITT)=ITT
126 : JPERM(ITT)=ITT
127 : *!!!!!!!!!!!!!!
128 : LWORK(ITT)=.FALSE.
129 : ENDDO
130 : *
131 : * On boucle sur les ddl
132 : *
133 : DO 1 ITT=1,NTT
134 : * Calcul de la norme de la ittème ligne
135 : ISTRT=IA(ITT)
136 : ISTOP=IA(ITT+1)-1
137 : TNORM=ZERO
138 : DO 11 IJA=ISTRT,ISTOP
139 : TNORM=TNORM+ABS(A(IJA))
140 : 11 CONTINUE
141 : TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
142 : IF (TNORM.LE.XPETIT) THEN
143 : WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
144 : GOTO 9999
145 : ENDIF
146 : * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
147 : * parties triangulaires inférieures et supérieurs.
148 : *
149 : * On peut le programmer autrement compte tenu du fait que les
150 : * colonnes sont ordonnées
151 : *
152 : LENU=1
153 : LENL=0
154 : IWORK(ITT)=ITT
155 : RWORK(ITT)=ZERO
156 : JWORK(ITT)=ITT
157 : LWORK(ITT)=.TRUE.
158 :
159 : DO 12 IJA=ISTRT,ISTOP
160 : *! COL=JA(IJA)
161 : COL=JPERM(JA(IJA))
162 : VAL=A(IJA)
163 : IF (COL.LT.ITT) THEN
164 : LENL=LENL+1
165 : IWORK(LENL)=COL
166 : RWORK(LENL)=VAL
167 : *!!!!!!!!!!!!!!
168 : LWORK(LENL)=.TRUE.
169 : JWORK(COL) =LENL
170 : ELSEIF (COL.EQ.ITT) THEN
171 : RWORK(ITT)=VAL
172 : ELSE
173 : LENU=LENU+1
174 : JPOS=ITT+LENU-1
175 : IWORK(JPOS)=COL
176 : RWORK(JPOS)=VAL
177 : *!!!!!!!!!!!!!!
178 : LWORK(JPOS)=.TRUE.
179 : JWORK(COL) =JPOS
180 : ENDIF
181 : 12 CONTINUE
182 : * On a enlevé les stratégies bizarres
183 : XNTLIG=DBLE(NNZA)/DBLE(NTT)
184 : *!! ILFIL=INT(XNTLIG*XLFIL/2.D0)
185 : ILFIL=INT(XNTLIG*(XLFIL-1.D0)/2.D0)
186 : ILFILL=ILFIL
187 : ILFILU=ILFIL
188 : * Elimination avec les lignes précédentes
189 : JJL=0
190 : LEN=0
191 : 13 CONTINUE
192 : JJL=JJL+1
193 : IF (JJL.LE.LENL) THEN
194 : c-----------------------------------------------------------------------
195 : c in order to do the elimination in the correct order we must select
196 : c the smallest column index among jw(k), k=jj+1, ..., lenl.
197 : c-----------------------------------------------------------------------
198 : JROW=IWORK(JJL)
199 : KJL=JJL
200 : c determine smallest column index
201 : DO 132 ILENL=JJL+1,LENL
202 : IF (IWORK(ILENL).LT.JROW) THEN
203 : JROW=IWORK(ILENL)
204 : KJL=ILENL
205 : ENDIF
206 : 132 CONTINUE
207 : IF (KJL.NE.JJL) THEN
208 : ITMP=IWORK(JJL)
209 : IWORK(JJL)=IWORK(KJL)
210 : IWORK(KJL)=ITMP
211 : JWORK(JROW)=JJL
212 : JWORK(ITMP)=KJL
213 : RTMP=RWORK(JJL)
214 : RWORK(JJL)=RWORK(KJL)
215 : RWORK(KJL)=RTMP
216 : *!!!!!!!!!!!!!!
217 : LTMP=LWORK(JJL)
218 : LWORK(JJL)=LWORK(KJL)
219 : LWORK(KJL)=LTMP
220 : ENDIF
221 : c zero out element in row by setting jw(n+jrow) to zero.
222 : JWORK(JROW)=0
223 : c get the multiplier for row to be eliminated (jrow).
224 : FACT=RWORK(JJL)*ALU(JROW)
225 : LFACT=LWORK(JJL)
226 : *!!!!!!!!!!!!!!
227 : IF (ABS(FACT).LE.XDTOL.AND.(.NOT.LFACT)) THEN
228 : * IF (ABS(FACT).LE.XDTOL) THEN
229 : GOTO 13
230 : ENDIF
231 : c
232 : c combine current row and row jrow
233 : c
234 : DO 134 KLU=JU(JROW),JLU(JROW+1)-1
235 : TERMLU=FACT*ALU(KLU)
236 : *! COLLU =JLU(KLU)
237 : COLLU =JPERM(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+1) THEN
245 : WRITE(IOIMP,*) 'pb. matrice U'
246 : GOTO 9999
247 : ENDIF
248 : IWORK(ITT+LENU-1)=COLLU
249 : JWORK(COLLU)=ITT+LENU-1
250 : RWORK(ITT+LENU-1)=-TERMLU
251 : LWORK(ITT+LENU-1)=.FALSE.
252 : ELSE
253 : c this is not a fill-in element
254 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
255 : ENDIF
256 : ELSE
257 : c dealing with lower part.
258 : IF (JPOS.EQ.0) THEN
259 : c this is a fill-in element
260 : LENL=LENL+1
261 : IF (LENL.GT.ITT-1) THEN
262 : WRITE(IOIMP,*) 'pb. matrice L'
263 : WRITE(IOIMP,*) 'LENL=',LENL
264 : WRITE(IOIMP,*) 'ITT=',ITT
265 : GOTO 9999
266 : ENDIF
267 : IWORK(LENL) =COLLU
268 : JWORK(COLLU)=LENL
269 : RWORK(LENL) =-TERMLU
270 : LWORK(LENL) =.FALSE.
271 : ELSE
272 : c this is not a fill-in element
273 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
274 : ENDIF
275 : ENDIF
276 : 134 CONTINUE
277 : c store this pivot element -- (from left to right -- no danger of
278 : c overlap with the working elements in L (pivots).
279 : LEN=LEN+1
280 : * ?????????
281 : LWORK(LEN)=LFACT
282 : RWORK(LEN)=FACT
283 : IWORK(LEN)=JROW
284 : GOTO 13
285 : ENDIF
286 : c reset double-pointer to zero (U-part including diagonal)
287 : DO 14 ILENU=1,LENU
288 : JWORK(IWORK(ITT+ILENU-1))=0
289 : 14 CONTINUE
290 : c update L-matrix
291 : LENL=LEN
292 : *
293 : * First store elements of the ilu(0) factorization
294 : *
295 : LEN2=0
296 : DO ILENL=1,LENL
297 : IF (LWORK(ILENL)) THEN
298 : LEN2=LEN2+1
299 : KWORK(LEN2)=IWORK(ILENL)
300 : QWORK(LEN2)=RWORK(ILENL)
301 : ENDIF
302 : ENDDO
303 : LENIL0=LEN2
304 : * then the fill-in elements
305 : DO ILENL=1,LENL
306 : IF (.NOT.LWORK(ILENL)) THEN
307 : LEN2=LEN2+1
308 : KWORK(LEN2)=IWORK(ILENL)
309 : QWORK(LEN2)=RWORK(ILENL)
310 : ENDIF
311 : ENDDO
312 : IF (LEN2.NE.LENL) THEN
313 : WRITE(IOIMP,*) 'Erreur de programmation 1'
314 : GOTO 9999
315 : ENDIF
316 : LENFIL=LENL-LENIL0
317 : * sort the fill-in elements by quick-split
318 : LEN=MIN(LENFIL,ILFILL)
319 : c sort by quick-split
319 : CALL QSPLIT(QWORK(LENIL0+1),KWORK(LENIL0+1),LENFIL,LEN,
320 : $ IMPR,IRET)
321 : c store L-part -- in original coordinates
322 : DO 15 ILEN=1,LENIL0+LEN
323 : ALU(JU0)=QWORK(ILEN)
324 : JLU(JU0)=IPERM(KWORK(ILEN))
325 : JU0=JU0+1
326 : 15 CONTINUE
327 : *
328 : c save pointer to beginning of row ii of U
329 : JU(ITT)=JU0
330 : c update U-matrix -- first apply dropping strategy
331 : LEN=0
332 : DO 16 ILEN=1,LENU-1
333 : IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
334 : LEN=LEN+1
335 : RWORK(ITT+LEN)=RWORK(ITT+ILEN)
336 : IWORK(ITT+LEN)=IWORK(ITT+ILEN)
337 : *!!!!!!!!!!!!!!
338 : LWORK(ITT+LEN)=LWORK(ITT+ILEN)
339 : ENDIF
340 : 16 CONTINUE
341 : LENU=LEN+1
342 : *
343 : * First store elements of the ilu(0) factorization
344 : *
345 : LEN2=ITT-1
346 : DO ILENU=ITT,ITT+LENU-1
347 : IF (LWORK(ILENU)) THEN
348 : LEN2=LEN2+1
349 : KWORK(LEN2)=IWORK(ILENU)
350 : QWORK(LEN2)=RWORK(ILENU)
351 : ENDIF
352 : ENDDO
353 : LENIL0=LEN2-ITT+1
354 : * then the fill-in elements
355 : DO ILENU=ITT,ITT+LENU-1
356 : IF (.NOT.LWORK(ILENU)) THEN
357 : LEN2=LEN2+1
358 : KWORK(LEN2)=IWORK(ILENU)
359 : QWORK(LEN2)=RWORK(ILENU)
360 : ENDIF
361 : ENDDO
362 : IF (LEN2-ITT+1.NE.LENU) THEN
363 : WRITE(IOIMP,*) 'Erreur de programmation 2'
364 : GOTO 9999
365 : ENDIF
366 : LENFIL=LENU-LENIL0
367 : * sort the fill-in elements by quicksplit
368 : *!!!!!!!!
369 : LEN=MIN(LENFIL,ILFILU)
370 : *!! LEN=MIN(LENU,ILFILU+1)
371 : *!!! LEN=MIN(LENU,ILFILU)
371 : CALL QSPLIT(QWORK(ITT+LENIL0),KWORK(ITT+LENIL0),LENFIL,LEN,
372 : $ IMPR,IRET)
373 : C WRITE(IOIMP,*) 'NTT=',NTT
374 : C WRITE(IOIMP,*) 'ITT=',ITT
375 : C WRITE(IOIMP,*) 'ILFILU=',ILFILU
376 : C WRITE(IOIMP,*) 'LEN=',LEN
377 : C WRITE(IOIMP,*) 'LENU=',LENU
378 : C DO JTT=ITT,NTT
379 : C WRITE(IOIMP,*) 'JTT=',JTT
380 : C WRITE(IOIMP,*) 'IWORK(JTT)=',IWORK(JTT)
381 : C WRITE(IOIMP,*) 'RWORK(JTT)=',RWORK(JTT)
382 : C WRITE(IOIMP,*) 'LWORK(JTT)=',LWORK(JTT)
383 : C ENDDO
384 : C DO JTT=ITT,NTT
385 : C WRITE(IOIMP,*) 'JTT=',JTT
386 : C WRITE(IOIMP,*) 'KWORK(JTT)=',KWORK(JTT)
387 : C WRITE(IOIMP,*) 'QWORK(JTT)=',QWORK(JTT)
388 : C ENDDO
389 : *
390 : * Determine next pivot
391 : *
392 : IMAX=ITT
393 : XMAX=ABS(QWORK(IMAX))
394 : IF (XSPIV.GT.ZERO) THEN
395 : XMAX=XMAX/XSPIV
396 : LPIV=.TRUE.
397 : * On pivote même si xspiv.EQ.zero si la diagonale est petite
398 : ELSEIF (XSPIV.EQ.ZERO.AND.XMAX.LE.PETIT*TNORM) THEN
399 : LPIV=.TRUE.
400 : ELSE
401 : LPIV=.FALSE.
402 : ENDIF
403 : *
404 : * WRITE(IOIMP,*) 'LPIV=',LPIV
405 : IF (LPIV) THEN
406 : DO ILEN=ITT+1,ITT+LENIL0+LEN-1
407 : XCOU=ABS(QWORK(ILEN))
408 : IF (XCOU.GT.XMAX) THEN
409 : IMAX=ILEN
410 : XMAX=XCOU
411 : ENDIF
412 : ENDDO
413 : * On ne cherche le max que dans le ilu0 sauf si on ne trouve rien
414 : C LFOUND=.FALSE.
415 : C DO ILEN=ITT+1,ITT+LENIL0-1
416 : C XCOU=ABS(QWORK(ILEN))
417 : C IF (XCOU.GT.XMAX) THEN
418 : C IMAX=ILEN
419 : C XMAX=XCOU
420 : C LFOUND=.TRUE.
421 : C ENDIF
422 : C ENDDO
423 : C IF (.NOT.LFOUND) THEN
424 : C DO ILEN=ITT+LENIL0,ITT+LENIL0+LEN-1
425 : C XCOU=ABS(QWORK(ILEN))
426 : C IF (XCOU.GT.XMAX) THEN
427 : C IMAX=ILEN
428 : C XMAX=XCOU
429 : C ENDIF
430 : C ENDDO
431 : C ENDIF
432 : *
433 : * Echange des valeurs
434 : *
435 : XTMP=QWORK(ITT)
436 : QWORK(ITT)=QWORK(IMAX)
437 : QWORK(IMAX)=XTMP
438 : *
439 : * Update de la numérotation ancienne
440 : *
441 : JMAX=KWORK(IMAX)
442 : ITMP=IPERM(ITT)
443 : IPERM(ITT)=IPERM(JMAX)
444 : IPERM(JMAX)=ITMP
445 : *
446 : * Update de la numérotation nouvelle
447 : *
448 : JPERM(IPERM(ITT))=ITT
449 : JPERM(IPERM(JMAX))=JMAX
450 : ENDIF
451 : *
452 : c copy U-Part in original coordinates
453 : *
454 : DO 17 ILEN=ITT+1,ITT+LENIL0+LEN-1
455 : JLU(JU0)=IPERM(KWORK(ILEN))
456 : ALU(JU0)=QWORK(ILEN)
457 : JU0=JU0+1
458 : 17 CONTINUE
459 : c store inverse of diagonal element of u
460 : VALPIV=QWORK(ITT)
461 : IF (VALPIV.EQ.ZERO) THEN
462 : WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
463 : $ 'le préconditionnement ILUTPGO2 est impossible.'
464 : IRET=-2
465 : GOTO 9999
466 : ELSE
467 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
468 : NBPIVP=NBPIVP+1
469 : ENDIF
470 : IF (VALPIV.LT.ZERO) THEN
471 : NBPIVI=NBPIVI+1
472 : ENDIF
473 : ENDIF
474 : ALU(ITT)=ONE/VALPIV
475 : c update pointer to beginning of next row of U.
476 : JLU(ITT+1)=JU0
477 : 1 CONTINUE
478 : NNZLU=JLU(NTT+1)-1
479 : C
480 : C Permutation des colonnes de LU
481 : C
482 : DO KLU=JLU(1),JLU(NTT+1)-1
483 : JLU(KLU)=JPERM(JLU(KLU))
484 : ENDDO
485 : C
486 : C Permutation des colonnes de A
487 : C
488 : DO IJA=IA(1),IA(NTT+1)-1
489 : JA(IJA)=JPERM(JA(IJA))
490 : ENDDO
491 : *
492 : * Warning(s)
493 : *
494 : IF (IMPR.GT.0) THEN
495 : IF (NBPIVP.GT.0) THEN
496 : WRITE(IOIMP,*) 'WARNING !',
497 : $ ' miltp3 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
498 : ENDIF
499 : *! IF (NBPIVI.GT.0) THEN
500 : *! WRITE(IOIMP,*) 'WARNING !',
501 : *! $ ' miltp3 :',NBPIVI,' pivot(s) < 0'
502 : *! ENDIF
503 : ENDIF
504 : *
505 : * Normal termination
506 : *
507 : IRET=0
508 : RETURN
509 : *
510 : * Error handling
511 : *
512 : 9999 CONTINUE
513 : IRET=1
514 : WRITE(IOIMP,*) 'An error was detected in miltp3.eso'
515 : RETURN
516 : *
517 : * End of MILTP3
518 : *
519 : END
520 :
521 :
522 :
523 :
524 :
525 :
© Cast3M 2003 - All rights reserved.
Disclaimer