1 : C MEILU0 SOURCE GOUNAND 05/02/16 21:16:06 5029
2 : SUBROUTINE MEILU0(NTT,NNZ,A,JA,IA,
3 : $ ALU,JLU,JU,
4 : $ IWORK,
5 : $ IMPR,IRET)
6 : IMPLICIT REAL*8 (A-H,O-Z)
7 : IMPLICIT INTEGER (I-N)
8 : C***********************************************************************
9 : C NOM : MEILU0
10 : C DESCRIPTION :
11 : C Calcul du préconditionneur ILU(0) d'une matrice Morse.
12 : C ILU(0) : Incomplete LU factorization of level 0
13 : C appelée aussi Choleski ou Crout incomplet
14 : C La matrice Morse est au format CSR (Column Sparse Row)
15 : C Le préconditionneur est au format MSR (Modified Sparse Row,
16 : C stockage de l'inverse de la diagonale) (ALU, JLU).
17 : C JU est un tableau supplémentaire pour le repérage des diagonales.
18 : C (voir la doc. de Sparskit version 2)
19 : C
20 : C LANGAGE : FORTRAN 77 (sauf E/S)
21 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
22 : C mél : gounand@semt2.smts.cea.fr
23 : C REFERENCE (bibtex-like) :
24 : C @BOOK{templates,
25 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
26 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
27 : C H. Van der Vorst},
28 : C TITLE={Templates for the Solution of Linear Systems :
29 : C Building Blocks for Iterative Methods},
30 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
31 : C -> URL : http://www.netlib.org/templates/Templates.html
32 : C Sparskit : a basic tool kit for sparse matrix computations
33 : C Version 2 (Youcef Saad)
34 : C -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
35 : C
36 : C***********************************************************************
37 : C APPELES : -
38 : C APPELE PAR : PRILU0
39 : C***********************************************************************
40 : C ENTREES : NTT, NNZ, IA, JA, A
41 : C ENTREES/SORTIES : IWORK (tableau de travail)
42 : C SORTIES : JU, JLU, ALU
43 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
44 : C***********************************************************************
45 : C VERSION : v1.1, 21/03/00
46 : C HISTORIQUE : v1.1, 22/03/00, gounand
47 : C Amélioration de la détection des pivots nulles.
48 : C HISTORIQUE : v1, 20/12/99, création
49 : C HISTORIQUE :
50 : C HISTORIQUE :
51 : C***********************************************************************
52 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
53 : C en cas de modification de ce sous-programme afin de faciliter
54 : C la maintenance !
55 : C***********************************************************************
56 : -INC CCOPTIO
57 : -INC CCREEL
58 : * .. Scalar Arguments ..
59 : INTEGER NTT,NNZ
60 : INTEGER IMPR,IRET
61 : * ..
62 : * .. Array Arguments ..
63 : REAL*8 A(NNZ),ALU(NNZ+1)
64 : INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
65 : INTEGER IWORK(NTT)
66 : *
67 : * .. Variables locales
68 : * .. Parameters
69 : REAL*8 ZERO ,ONE
70 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
71 : * ..
72 : C Nombre de pivots petits
73 : INTEGER NBPIVP
74 : C Nombre de pivots inférieurs à 0
75 : INTEGER NBPIVI
76 : INTEGER ITT,JNZ,KNZ,JCOL,JROW
77 : INTEGER JSTRT,JMID,JSTOP,JU0,JW
78 : INTEGER ISTRT,ISTOP
79 : REAL*8 TL,VALPIV,TNORM,PETIT
80 :
81 : *
82 : * Executable statements
83 : *
84 : PETIT=XZPREC**(0.75D0)
85 : NBPIVP=0
86 : NBPIVI=0
87 : JU0=NTT+2
88 : JLU(1)=JU0
89 : *
90 : * On boucle sur les ddl
91 : *
92 : DO 1 ITT=1,NTT
93 : JSTRT=JU0
94 : *
95 : * Création de la ligne ITT pour L et U
96 : *
97 : ISTRT=IA(ITT)
98 : ISTOP=IA(ITT+1)-1
99 : TNORM=ZERO
100 : DO 12 JNZ=ISTRT,ISTOP
101 : *
102 : * Copie de la ligne ITT du format CSR au format MSR
103 : * et calcul de la norme de la ligne
104 : *
105 : JCOL=JA(JNZ)
106 : TNORM=TNORM+ABS(A(JNZ))
107 : IF (JCOL.EQ.ITT) THEN
108 : ALU(ITT) =A(JNZ)
109 : IWORK(JCOL)=ITT
110 : JU(ITT) =JU0
111 : ELSE
112 : ALU(JU0) =A(JNZ)
113 : JLU(JU0) =JA(JNZ)
114 : IWORK(JCOL)=JU0
115 : JU0 = JU0+1
116 : ENDIF
117 : 12 CONTINUE
118 : TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
119 : IF (TNORM.LE.XPETIT) THEN
120 : WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
121 : GOTO 9999
122 : ENDIF
123 : JLU(ITT+1)=JU0
124 : JSTOP=JU0-1
125 : JMID=JU(ITT)-1
126 : DO 14 JNZ=JSTRT,JMID
127 : JROW=JLU(JNZ)
128 : TL=ALU(JNZ)*ALU(JROW)
129 : ALU(JNZ)=TL
130 : *
131 : * On effectue la combinaison linéaire
132 : *
133 : DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
134 : JW=IWORK(JLU(KNZ))
135 : IF (JW.NE.0) THEN
136 : ALU(JW)=ALU(JW)-TL*ALU(KNZ)
137 : ENDIF
138 : 142 CONTINUE
139 : 14 CONTINUE
140 : *
141 : * On s'occupe de l'élément diagonal
142 : *
143 : VALPIV=ALU(ITT)
144 : IF (VALPIV.EQ.ZERO) THEN
145 : WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
146 : $ 'le préconditionnement ILU(0) est impossible.'
147 : IRET=-2
148 : GOTO 9999
149 : ELSE
150 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
151 : NBPIVP=NBPIVP+1
152 : ENDIF
153 : IF (VALPIV.LT.ZERO) THEN
154 : NBPIVI=NBPIVI+1
155 : ENDIF
156 : ENDIF
157 : ALU(ITT)=ONE/ALU(ITT)
158 : *
159 : * On remet le tableau de travail à zéro
160 : *
161 : IWORK(ITT)=0
162 : DO 16 JNZ=JSTRT,JSTOP
163 : IWORK(JLU(JNZ))=0
164 : 16 CONTINUE
165 : 1 CONTINUE
166 : *
167 : * Warning(s)
168 : *
169 : IF (IMPR.GT.0) THEN
170 : IF (NBPIVP.GT.0) THEN
171 : WRITE(IOIMP,*) 'WARNING !',
172 : $ ' meilu0 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
173 : ENDIF
174 : *! IF (NBPIVI.GT.0) THEN
175 : *! WRITE(IOIMP,*) 'WARNING !',
176 : *! $ ' meilu0 :',NBPIVI,' pivot(s) < 0'
177 : *! ENDIF
178 : ENDIF
179 : *
180 : * Normal termination
181 : *
182 : IRET=0
183 : RETURN
184 : *
185 : * Error handling
186 : *
187 : 9999 CONTINUE
188 : IRET=1
189 : WRITE(IOIMP,*) 'An error was detected in meilu0.eso'
190 : RETURN
191 : *
192 : * End of MEILU0
193 : *
194 : END
195 :
196 :
197 :
198 :
199 :
200 :
© Cast3M 2003 - All rights reserved.
Disclaimer