1 : C MEMILU SOURCE GOUNAND 06/03/06 21:18:28 5319
2 : SUBROUTINE MEMILU(NTT,NNZ,A,JA,IA,
3 : $ ALU,JLU,JU,
4 : $ IWORK,
5 : $ RXMILU,
6 : $ IMPR,IRET)
7 : IMPLICIT INTEGER(I-N)
8 : IMPLICIT REAL*8 (A-H,O-Z)
9 : C***********************************************************************
10 : C NOM : MEMILU
11 : C DESCRIPTION :
12 : C Calcul du préconditionneur MILU(0) d'une matrice Morse
13 : C MILU(0) : Modified Incomplete LU factorization of level 0
14 : C appelée aussi Choleski ou Crout incomplet modifié !
15 : C Ici, c'est meme un MILU(0) relaxé (avec RXMILU)...
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, RXMILU
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 RXMILU : type REAL*8.
47 : C facteur de relaxation (normalement compris
48 : C entre 0.D0 et 1.D0)
49 : C 0.D0 : on se ramène à Choleski incomplet
50 : C 1.D0 : on se ramène à MILU(0) non relaxé.
51 : C***********************************************************************
52 : C VERSION : v1.1, 21/03/00
53 : C HISTORIQUE : v1.1, 22/03/00, gounand
54 : C Amélioration de la détection des pivots nulles.
55 : C HISTORIQUE : v1, 20/12/99, 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,NNZ
67 : REAL*8 RXMILU
68 : INTEGER IMPR,IRET
69 : * ..
70 : * .. Array Arguments ..
71 : REAL*8 A(NNZ),ALU(NNZ+1)
72 : INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
73 : INTEGER IWORK(NNZ)
74 : *
75 : * .. Variables locales
76 : * .. Parameters
77 : REAL*8 ZERO ,ONE
78 : PARAMETER (ZERO=0.0D0,ONE=1.0D0)
79 : * ..
80 : C Nombre de pivots petits
81 : INTEGER NBPIVP
82 : C Nombre de pivots inférieurs à 0
83 : INTEGER NBPIVI
84 : INTEGER ITT,JNZ,KNZ,JCOL,JROW
85 : INTEGER JSTRT,JMID,JSTOP,JU0,JW
86 : INTEGER ISTRT,ISTOP
87 : REAL*8 TL,VALPIV,TNORM,PETIT
88 : REAL*8 FILVAL
89 : *
90 : * Executable statements
91 : *
92 : PETIT=XZPREC**(0.75D0)
93 : NBPIVP=0
94 : NBPIVI=0
95 : JU0=NTT+2
96 : JLU(1)=JU0
97 : *
98 : * On boucle sur les ddl
99 : *
100 : DO 1 ITT=1,NTT
101 : JSTRT=JU0
102 : *
103 : * Création de la ligne ITT pour L et U
104 : *
105 : ISTRT=IA(ITT)
106 : ISTOP=IA(ITT+1)-1
107 : TNORM=ZERO
108 : DO 12 JNZ=ISTRT,ISTOP
109 : *
110 : * Copie de la ligne ITT du format CSR au format MSR
111 : *
112 : JCOL=JA(JNZ)
113 : TNORM=TNORM+ABS(A(JNZ))
114 : IF (JCOL.EQ.ITT) THEN
115 : ALU(ITT) =A(JNZ)
116 : IWORK(JCOL)=ITT
117 : JU(ITT) =JU0
118 : ELSE
119 : ALU(JU0) =A(JNZ)
120 : JLU(JU0) =JA(JNZ)
121 : IWORK(JCOL)=JU0
122 : JU0 = JU0+1
123 : ENDIF
124 : 12 CONTINUE
125 : TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
126 : IF (TNORM.LE.XPETIT) THEN
127 : WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
128 : GOTO 9999
129 : ENDIF
130 : JLU(ITT+1)=JU0
131 : JSTOP=JU0-1
132 : JMID=JU(ITT)-1
133 : FILVAL=ZERO
134 : DO 14 JNZ=JSTRT,JMID
135 : JROW=JLU(JNZ)
136 : TL=ALU(JNZ)*ALU(JROW)
137 : ALU(JNZ)=TL
138 : *
139 : * On effectue la combinaison linéaire
140 : *
141 : DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
142 : JW=IWORK(JLU(KNZ))
143 : IF (JW.NE.0) THEN
144 : ALU(JW)=ALU(JW)-TL*ALU(KNZ)
145 : ELSE
146 : FILVAL=FILVAL+TL*ALU(KNZ)
147 : ENDIF
148 : 142 CONTINUE
149 : 14 CONTINUE
150 : *
151 : * On s'occupe de l'élément diagonal
152 : *
153 : ALU(ITT)=ALU(ITT)-(FILVAL*RXMILU)
154 : VALPIV=ALU(ITT)
155 : IF (VALPIV.EQ.ZERO) THEN
156 : WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
157 : $ 'le préconditionnement MILU(0) est impossible.'
158 : IRET=-2
159 : GOTO 9999
160 : ELSE
161 : IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
162 : NBPIVP=NBPIVP+1
163 : ENDIF
164 : IF (VALPIV.LT.ZERO) THEN
165 : NBPIVI=NBPIVI+1
166 : ENDIF
167 : ENDIF
168 : ALU(ITT)=ONE/ALU(ITT)
169 : *
170 : * On remet le tableau de travail à zéro
171 : *
172 : IWORK(ITT)=0
173 : DO 16 JNZ=JSTRT,JSTOP
174 : IWORK(JLU(JNZ))=0
175 : 16 CONTINUE
176 : 1 CONTINUE
177 : *
178 : * Warning(s)
179 : *
180 : IF (IMPR.GT.0) THEN
181 : IF (NBPIVP.GT.0) THEN
182 : WRITE(IOIMP,*) 'WARNING !',
183 : $ ' memilu :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
184 : ENDIF
185 : * IF (NBPIVI.GT.0) THEN
186 : * WRITE(IOIMP,*) 'WARNING !',
187 : * $ ' memilu :',NBPIVI,' pivot(s) < 0'
188 : * ENDIF
189 : ENDIF
190 : *
191 : * Normal termination
192 : *
193 : IRET=0
194 : RETURN
195 : *
196 : * Error handling
197 : *
198 : 9999 CONTINUE
199 : WRITE(IOIMP,*) 'An error was detected in memilu.eso'
200 : RETURN
201 : *
202 : * End of MEMILU
203 : *
204 : END
205 :
206 :
207 :
208 :
209 :
210 :
211 :
© Cast3M 2003 - All rights reserved.
Disclaimer