1 : C MEILUP SOURCE GOUNAND 11/07/21 21:15:25 7046
2 : SUBROUTINE MEILUP(NTT,NNZ,A,JA,IA,
3 : $ ALU,JLU,JU,
4 : $ IWORK,
5 : $ RXILUP,IFILTR,
6 : $ IMPR,IRET)
7 : IMPLICIT REAL*8 (A-H,O-Z)
8 : IMPLICIT INTEGER (I-N)
9 : C***********************************************************************
10 : C NOM : MEILUP
11 : C DESCRIPTION :
12 : C Calcul du préconditionneur ILU(0) d'une matrice Morse.
13 : C modifié PV
14 : C ILU(0) : Incomplete LU factorization of level 0
15 : C appelée aussi Choleski ou Crout incomplet
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
38 : C***********************************************************************
39 : C APPELES : -
40 : C APPELE PAR : PRILU0
41 : C***********************************************************************
42 : C ENTREES : NTT, NNZ, IA, JA, A
43 : C ENTREES/SORTIES : IWORK (tableau de travail)
44 : C SORTIES : JU, JLU, ALU
45 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
46 : C***********************************************************************
47 : C VERSION : v1.1, 21/03/00
48 : C HISTORIQUE : v1.1, 22/03/00, gounand
49 : C Amélioration de la détection des pivots nulles.
50 : C HISTORIQUE : v1, 20/12/99, création
51 : C HISTORIQUE :
52 : C HISTORIQUE :
53 : C***********************************************************************
54 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
55 : C en cas de modification de ce sous-programme afin de faciliter
56 : C la maintenance !
57 : C***********************************************************************
58 : -INC CCOPTIO
59 : -INC CCREEL
60 : * .. Scalar Arguments ..
61 : INTEGER NTT,NNZ
62 : INTEGER IMPR,IRET
63 : * ..
64 : * .. Array Arguments ..
65 : REAL*8 A(NNZ),ALU(NNZ+1)
66 : INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
67 : INTEGER IWORK(NTT)
68 : *
69 : * .. Variables locales
70 : * .. Parameters
71 : REAL*8 ZERO ,ONE
72 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
73 : * ..
74 : C Nombre de pivots petits
75 : INTEGER NBPIVP
76 : C Nombre de pivots inférieurs à 0
77 : INTEGER NBPIVI
78 : INTEGER ITT,JNZ,KNZ,JCOL,JROW
79 : INTEGER JSTRT,JMID,JSTOP,JU0,JW
80 : INTEGER ISTRT,ISTOP
81 : REAL*8 TL,VALPIV,TNORM,PETIT,XOPIV,XNPIV
82 :
83 : *
84 : * Executable statements
85 : *
86 : PETIT=XZPREC**(0.75D0)
87 : NBPIVP=0
88 : NBPIVI=0
89 : JU0=NTT+2
90 : JLU(1)=JU0
91 : *
92 : * On boucle sur les ddl
93 : *
94 : DO 1 ITT=1,NTT
95 : JSTRT=JU0
96 : *
97 : * Création de la ligne ITT pour L et U
98 : *
99 : ISTRT=IA(ITT)
100 : ISTOP=IA(ITT+1)-1
101 : TNORM=ZERO
102 : XOPIV=ZERO
103 : DO 12 JNZ=ISTRT,ISTOP
104 : *
105 : * Copie de la ligne ITT du format CSR au format MSR
106 : * et calcul de la norme de la ligne
107 : *
108 : JCOL=JA(JNZ)
109 : TNORM=TNORM+ABS(A(JNZ))
110 : IF (JCOL.EQ.ITT) THEN
111 : XOPIV =A(JNZ)
112 : ALU(ITT) =A(JNZ)
113 : IWORK(JCOL)=ITT
114 : JU(ITT) =JU0
115 : ELSE
116 : ALU(JU0) =A(JNZ)
117 : JLU(JU0) =JA(JNZ)
118 : IWORK(JCOL)=JU0
119 : JU0 = JU0+1
120 : ENDIF
121 : 12 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 : JLU(ITT+1)=JU0
128 : JSTOP=JU0-1
129 : JMID=JU(ITT)-1
130 : DO 14 JNZ=JSTRT,JMID
131 : JROW=JLU(JNZ)
132 : TL=ALU(JNZ)*ALU(JROW)
133 : ALU(JNZ)=TL
134 : *
135 : * On effectue la combinaison linéaire
136 : *
137 : DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
138 : JW=IWORK(JLU(KNZ))
139 : IF (JW.NE.0) THEN
140 : ALU(JW)=ALU(JW)-TL*ALU(KNZ)
141 : ENDIF
142 : 142 CONTINUE
143 : 14 CONTINUE
144 : *
145 : * On s'occupe de l'élément diagonal
146 : *
147 : *ILU(0) VALPIV=ALU(ITT)
148 : IF (IFILTR.EQ.1) THEN
149 : *New deb
150 : XNPIV=ALU(ITT)
151 : VALPIV=XNPIV
152 : IF (ABS(XOPIV).GT.XPETIT) THEN
153 : IF (XNPIV.LT.(ABS(XOPIV)*RXILUP)) THEN
154 : VALPIV=XOPIV
155 : ENDIF
156 : ENDIF
157 : ENDIF
158 : *New fin
159 : IF (IFILTR.EQ.0) THEN
160 : VALPIV=XOPIV
161 : ENDIF
162 : *
163 : IF (VALPIV.EQ.ZERO) THEN
164 : WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
165 : $ 'le préconditionnement ILU(0)-PV est impossible.'
166 : IRET=-2
167 : GOTO 9999
168 : ELSE
169 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
170 : NBPIVP=NBPIVP+1
171 : ENDIF
172 : IF (VALPIV.LT.ZERO) THEN
173 : NBPIVI=NBPIVI+1
174 : ENDIF
175 : ENDIF
176 : ALU(ITT)=ONE/VALPIV
177 : *
178 : * On remet le tableau de travail à zéro
179 : *
180 : IWORK(ITT)=0
181 : DO 16 JNZ=JSTRT,JSTOP
182 : IWORK(JLU(JNZ))=0
183 : 16 CONTINUE
184 : 1 CONTINUE
185 : *
186 : * Warning(s)
187 : *
188 : IF (IMPR.GT.0) THEN
189 : IF (NBPIVP.GT.0) THEN
190 : WRITE(IOIMP,*) 'WARNING !',
191 : $ ' meilup :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
192 : ENDIF
193 : *! IF (NBPIVI.GT.0) THEN
194 : *! WRITE(IOIMP,*) 'WARNING !',
195 : *! $ ' meilup :',NBPIVI,' pivot(s) < 0'
196 : *! ENDIF
197 : ENDIF
198 : *
199 : * Normal termination
200 : *
201 : IRET=0
202 : RETURN
203 : *
204 : * Error handling
205 : *
206 : 9999 CONTINUE
207 : IRET=1
208 : WRITE(IOIMP,*) 'An error was detected in meilup.eso'
209 : RETURN
210 : *
211 : * End of MEILUP
212 : *
213 : END
214 :
215 :
216 :
217 :
218 :
219 :
220 :
© Cast3M 2003 - All rights reserved.
Disclaimer