1 : C MEILTP SOURCE GOUNAND 05/06/06 21:15:36 5078
2 : SUBROUTINE MEILTP(NTT,NNZA,A,JA,IA,
3 : $ XLFIL,XDTOL,XSPIV,NNZMLU,
4 : $ ALU,JLU,JU,NNZLU,
5 : $ IWORK,JWORK,RWORK,
6 : $ IPERM,JPERM,
7 : $ IMPR,IRET)
8 : IMPLICIT REAL*8 (A-H,O-Z)
9 : IMPLICIT INTEGER (I-N)
10 : C***********************************************************************
11 : C NOM : MEILTP
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 a tenté de corriger les bugs de la version Saad
43 : C signalés par *!
44 : C
45 : C***********************************************************************
46 : C APPELES : QSPLIT
47 : C APPELE PAR : PRILUT
48 : C***********************************************************************
49 : C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
50 : C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
51 : C SORTIES : ALU,JLU,JU,NNZLU
52 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
53 : C***********************************************************************
54 : C VERSION : v1, 08/04/04, version initiale
55 : C HISTORIQUE : v1, 08/04/04, création
56 : C HISTORIQUE :
57 : C HISTORIQUE :
58 : C***********************************************************************
59 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
60 : C en cas de modification de ce sous-programme afin de faciliter
61 : C la maintenance !
62 : C***********************************************************************
63 : -INC CCOPTIO
64 : -INC CCREEL
65 : * .. Scalar Arguments ..
66 : INTEGER NTT,NNZA,NNZMLU,NNZLU
67 : INTEGER IMPR,IRET
68 : REAL*8 XDTOL,XLFIL
69 : * ..
70 : * .. Array Arguments ..
71 : REAL*8 A(NNZA),ALU(NNZMLU+1)
72 : INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
73 : INTEGER IWORK(NTT)
74 : INTEGER JWORK(NTT)
75 : REAL*8 RWORK(NTT)
76 : INTEGER IPERM(NTT)
77 : INTEGER JPERM(NTT)
78 : *
79 : * .. Variables locales
80 : * .. Parameters
81 : REAL*8 ZERO ,ONE
82 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
83 : * ..
84 : C Nombre de pivots petits
85 : INTEGER NBPIVP
86 : C Nombre de pivots inférieurs à 0
87 : INTEGER NBPIVI
88 : *
89 : INTEGER ILFIL,ILFILL,ILFILU
90 : INTEGER COL,COLLU
91 : INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
92 : INTEGER JJL,JPOS,JROW,KJL,KLU
93 : INTEGER LEN,LENL,LENU
94 : INTEGER JU0,ITMP
95 : REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
96 : REAL*8 XNTLIG
97 : * REAL*8 XPOND1,XPOND2
98 : *
99 : * Executable statements
100 : *
101 : IF (XLFIL.LT.0.D0) THEN
102 : WRITE(IOIMP,*) 'XLFIL.LT.0.D0 non acceptable'
103 : GOTO 9999
104 : ENDIF
105 : PETIT=XZPREC**(0.75D0)
106 : NBPIVP=0
107 : NBPIVI=0
108 : c-----------------------------------------------------------------------
109 : c initialize ju0 (points to next element to be added to alu,jlu)
110 : c and pointer array.
111 : c-----------------------------------------------------------------------
112 : JU0=NTT+2
113 : JLU(1)=JU0
114 : *
115 : * Initialize permutation arrays
116 : *
117 : DO ITT=1,NTT
118 : JWORK(ITT)=0
119 : IPERM(ITT)=ITT
120 : JPERM(ITT)=ITT
121 : ENDDO
122 : *
123 : * On boucle sur les ddl
124 : *
125 : DO 1 ITT=1,NTT
126 : * Calcul de la norme de la ittème ligne
127 : ISTRT=IA(ITT)
128 : ISTOP=IA(ITT+1)-1
129 : TNORM=ZERO
130 : DO 11 IJA=ISTRT,ISTOP
131 : TNORM=TNORM+ABS(A(IJA))
132 : 11 CONTINUE
133 : TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
134 : IF (TNORM.LE.XPETIT) THEN
135 : WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
136 : GOTO 9999
137 : ENDIF
138 : * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
139 : * parties triangulaires inférieures et supérieurs.
140 : *
141 : * On peut le programmer autrement compte tenu du fait que les
142 : * colonnes sont ordonnées
143 : *
144 : LENU=1
145 : LENL=0
146 : IWORK(ITT)=ITT
147 : RWORK(ITT)=ZERO
148 : JWORK(ITT)=ITT
149 :
150 : DO 12 IJA=ISTRT,ISTOP
151 : *! COL=JA(IJA)
152 : COL=JPERM(JA(IJA))
153 : VAL=A(IJA)
154 : IF (COL.LT.ITT) THEN
155 : LENL=LENL+1
156 : IWORK(LENL)=COL
157 : RWORK(LENL)=VAL
158 : JWORK(COL) =LENL
159 : ELSEIF (COL.EQ.ITT) THEN
160 : RWORK(ITT)=VAL
161 : ELSE
162 : LENU=LENU+1
163 : JPOS=ITT+LENU-1
164 : IWORK(JPOS)=COL
165 : RWORK(JPOS)=VAL
166 : JWORK(COL) =JPOS
167 : ENDIF
168 : 12 CONTINUE
169 : * On a enlevé les stratégies bizarres
170 : XNTLIG=DBLE(NNZA)/DBLE(NTT)
171 : ILFIL=INT((XNTLIG-1.D0)*XLFIL/2.D0)
172 : ILFILL=ILFIL
173 : ILFILU=ILFIL
174 : * Elimination avec les lignes précédentes
175 : JJL=0
176 : LEN=0
177 : 13 CONTINUE
178 : JJL=JJL+1
179 : IF (JJL.LE.LENL) THEN
180 : c-----------------------------------------------------------------------
181 : c in order to do the elimination in the correct order we must select
182 : c the smallest column index among jw(k), k=jj+1, ..., lenl.
183 : c-----------------------------------------------------------------------
184 : JROW=IWORK(JJL)
185 : KJL=JJL
186 : c determine smallest column index
187 : DO 132 ILENL=JJL+1,LENL
188 : IF (IWORK(ILENL).LT.JROW) THEN
189 : JROW=IWORK(ILENL)
190 : KJL=ILENL
191 : ENDIF
192 : 132 CONTINUE
193 : IF (KJL.NE.JJL) THEN
194 : ITMP=IWORK(JJL)
195 : IWORK(JJL)=IWORK(KJL)
196 : IWORK(KJL)=ITMP
197 : JWORK(JROW)=JJL
198 : JWORK(ITMP)=KJL
199 : RTMP=RWORK(JJL)
200 : RWORK(JJL)=RWORK(KJL)
201 : RWORK(KJL)=RTMP
202 : ENDIF
203 : c zero out element in row by setting jw(n+jrow) to zero.
204 : JWORK(JROW)=0
205 : c get the multiplier for row to be eliminated (jrow).
206 : FACT=RWORK(JJL)*ALU(JROW)
207 : IF (ABS(FACT).LE.XDTOL) THEN
208 : GOTO 13
209 : ENDIF
210 : c
211 : c combine current row and row jrow
212 : c
213 : DO 134 KLU=JU(JROW),JLU(JROW+1)-1
214 : TERMLU=FACT*ALU(KLU)
215 : *! COLLU =JLU(KLU)
216 : COLLU =JPERM(JLU(KLU))
217 : JPOS =JWORK(COLLU)
218 : IF (COLLU.GE.ITT) THEN
219 : c dealing with upper part (including diagonal)
220 : IF (JPOS.EQ.0) THEN
221 : c this is a fill-in element
222 : LENU=LENU+1
223 : IF (LENU.GT.NTT-ITT+1) THEN
224 : WRITE(IOIMP,*) 'pb. matrice U'
225 : GOTO 9999
226 : ENDIF
227 : IWORK(ITT+LENU-1)=COLLU
228 : JWORK(COLLU)=ITT+LENU-1
229 : RWORK(ITT+LENU-1)=-TERMLU
230 : ELSE
231 : c this is not a fill-in element
232 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
233 : ENDIF
234 : ELSE
235 : c dealing with lower part.
236 : IF (JPOS.EQ.0) THEN
237 : c this is a fill-in element
238 : LENL=LENL+1
239 : IF (LENL.GT.ITT-1) THEN
240 : WRITE(IOIMP,*) 'pb. matrice L'
241 : WRITE(IOIMP,*) 'LENL=',LENL
242 : WRITE(IOIMP,*) 'ITT=',ITT
243 : GOTO 9999
244 : ENDIF
245 : IWORK(LENL) =COLLU
246 : JWORK(COLLU)=LENL
247 : RWORK(LENL) =-TERMLU
248 : ELSE
249 : c this is not a fill-in element
250 : RWORK(JPOS)=RWORK(JPOS)-TERMLU
251 : ENDIF
252 : ENDIF
253 : 134 CONTINUE
254 : c store this pivot element -- (from left to right -- no danger of
255 : c overlap with the working elements in L (pivots).
256 : LEN=LEN+1
257 : RWORK(LEN)=FACT
258 : IWORK(LEN)=JROW
259 : GOTO 13
260 : ENDIF
261 : c reset double-pointer to zero (U-part including diagonal)
262 : DO 14 ILENU=1,LENU
263 : JWORK(IWORK(ITT+ILENU-1))=0
264 : 14 CONTINUE
265 : c update L-matrix
266 : LENL=LEN
267 : LEN=MIN(LENL,ILFILL)
268 : c sort by quick-split
268 : CALL QSPLIT(RWORK,IWORK,LENL,LEN,
269 : $ IMPR,IRET)
270 : c store L-part -- in original coordinates
271 : DO 15 ILEN=1,LEN
272 : *! JLU(JU0)=IWORK(ILEN)
273 : ALU(JU0)=RWORK(ILEN)
274 : JLU(JU0)=IPERM(IWORK(ILEN))
275 : JU0=JU0+1
276 : 15 CONTINUE
277 : c save pointer to beginning of row ii of U
278 : JU(ITT)=JU0
279 : c update U-matrix -- first apply dropping strategy
280 : LEN=0
281 : DO 16 ILEN=1,LENU-1
282 : IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
283 : LEN=LEN+1
284 : RWORK(ITT+LEN)=RWORK(ITT+ILEN)
285 : IWORK(ITT+LEN)=IWORK(ITT+ILEN)
286 : ENDIF
287 : 16 CONTINUE
288 : LENU=LEN+1
289 : *!!!!!!!!
290 : LEN=MIN(LENU,ILFILU+1)
291 : *!!! LEN=MIN(LENU,ILFILU)
292 : c sort by quick-split
292 : CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU-1,LEN-1,
293 : $ IMPR,IRET)
294 : *!
295 : *! Determine next pivot
296 : *!
297 : IMAX=ITT
298 : XMAX=ABS(RWORK(IMAX))
299 : XMAX0=XMAX
300 : * Je ne comprends pas bien cette ligne, je ne la mets pas
301 : * ICUT=ITT-1+IRPIV-MOD(ITT-1,IRPIV)
302 : * ICUT=ITT+IRPIV
303 : DO ILEN=ITT+1,ITT+LEN-1
304 : *!!! DO ILEN=ITT+1,ITT+LEN
305 : XCOU=ABS(RWORK(ILEN))
306 : IF (XCOU.GT.XMAX.AND.(XCOU*XSPIV).GT.XMAX0) THEN
307 : IMAX=ILEN
308 : XMAX=XCOU
309 : ENDIF
310 : ENDDO
311 : *
312 : * Echange des valeurs
313 : *
314 : XTMP=RWORK(ITT)
315 : RWORK(ITT)=RWORK(IMAX)
316 : RWORK(IMAX)=XTMP
317 : *
318 : * Update de la numérotation ancienne
319 : *
320 : JMAX=IWORK(IMAX)
321 : ITMP=IPERM(ITT)
322 : IPERM(ITT)=IPERM(JMAX)
323 : IPERM(JMAX)=ITMP
324 : *
325 : * Update de la numérotation nouvelle
326 : *
327 : JPERM(IPERM(ITT))=ITT
328 : JPERM(IPERM(JMAX))=JMAX
329 : *
330 : c copy U-Part in original coordinates
331 : *
332 : DO 17 ILEN=ITT+1,ITT+LEN-1
333 : *!!! DO 17 ILEN=ITT+1,ITT+LEN
334 : JLU(JU0)=IPERM(IWORK(ILEN))
335 : ALU(JU0)=RWORK(ILEN)
336 : JU0=JU0+1
337 : 17 CONTINUE
338 : c store inverse of diagonal element of u
339 : VALPIV=RWORK(ITT)
340 : IF (VALPIV.EQ.ZERO) THEN
341 : WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
342 : $ 'le préconditionnement ILUTP est impossible.'
343 : IRET=-2
344 : GOTO 9999
345 : ELSE
346 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
347 : NBPIVP=NBPIVP+1
348 : ENDIF
349 : IF (VALPIV.LT.ZERO) THEN
350 : NBPIVI=NBPIVI+1
351 : ENDIF
352 : ENDIF
353 : ALU(ITT)=ONE/RWORK(ITT)
354 : c update pointer to beginning of next row of U.
355 : JLU(ITT+1)=JU0
356 : 1 CONTINUE
357 : NNZLU=JLU(NTT+1)-1
358 : C
359 : C Permutation des colonnes de LU
360 : C
361 : DO KLU=JLU(1),JLU(NTT+1)-1
362 : JLU(KLU)=JPERM(JLU(KLU))
363 : ENDDO
364 : C
365 : C Permutation des colonnes de A
366 : C
367 : DO IJA=IA(1),IA(NTT+1)-1
368 : JA(IJA)=JPERM(JA(IJA))
369 : ENDDO
370 : *
371 : * Warning(s)
372 : *
373 : IF (IMPR.GT.0) THEN
374 : IF (NBPIVP.GT.0) THEN
375 : WRITE(IOIMP,*) 'WARNING !',
376 : $ ' meiltp :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
377 : ENDIF
378 : *! IF (NBPIVI.GT.0) THEN
379 : *! WRITE(IOIMP,*) 'WARNING !',
380 : *! $ ' meiltp :',NBPIVI,' pivot(s) < 0'
381 : *! ENDIF
382 : ENDIF
383 : *
384 : * Normal termination
385 : *
386 : IRET=0
387 : RETURN
388 : *
389 : * Error handling
390 : *
391 : 9999 CONTINUE
392 : IRET=1
393 : WRITE(IOIMP,*) 'An error was detected in meiltp.eso'
394 : RETURN
395 : *
396 : * End of MEILTP
397 : *
398 : END
399 :
400 :
401 :
402 :
403 :
404 :
405 :
© Cast3M 2003 - All rights reserved.
Disclaimer