1 : C MDILUN SOURCE GOUNAND 06/03/06 21:18:10 5319
2 : SUBROUTINE MDILUN(N,NNZ,
3 : $ ROWPTR,COLIND,VAL,
4 : $ DINV,
5 : $ IRET)
6 : IMPLICIT INTEGER(I-N)
7 : IMPLICIT REAL*8 (A-H,O-Z)
8 : C***********************************************************************
9 : C NOM : MDILUN
10 : C DESCRIPTION :
11 : C Calcul du préconditionneur D-ILU d'une matrice Morse
12 : C NON symétrique.
13 : C D-ILU : Diagonal Incomplete LU factorization
14 : C
15 : C Ce sous-programme est appelé par prdilu.eso.
16 : C Il est écrit en Fortran car il est censé etre rapide
17 : C
18 : C LANGAGE : FORTRAN 77 + chouilla ESOPE (E/S)
19 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
20 : C mél : gounand@semt2.smts.cea.fr
21 : C REFERENCE (bibtex-like) :
22 : C @BOOK{templates,
23 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
24 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
25 : C H. Van der Vorst},
26 : C TITLE={Templates for the Solution of Linear Systems :
27 : C Building Blocks for Iterative Methods},
28 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
29 : C -> URL : http://www.netlib.org/templates/Templates.html
30 : C***********************************************************************
31 : C APPELES : -
32 : C***********************************************************************
33 : C ENTREES : N, NNZ,
34 : C ROWPTR, COLIND, VAL
35 : C ENTREES/SORTIES : -
36 : C SORTIES : DINVPIV, IRET
37 : C CODE RETOUR (IRET) : 0 si ok
38 : C <0 si problème (diagonale nulle...)
39 : C N : nombre de degrés de liberté
40 : C NNZ : nombre de valeurs non nulles de la matrice Morse
41 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
42 : C et valeurs de la matrice Morse
43 : C DINVPIV : vecteur contenant l'inverse de la factorisation
44 : C D-ILU de la matrice morse ci-dessus.
45 : C***********************************************************************
46 : C VERSION : v1, 01/04/98, version initiale
47 : C HISTORIQUE : v1, 01/04/98, création
48 : C HISTORIQUE :
49 : C HISTORIQUE :
50 : C***********************************************************************
51 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
52 : C en cas de modification de ce sous-programme afin de faciliter
53 : C la maintenance !
54 : C***********************************************************************
55 : -INC CCOPTIO
56 : * .. Scalar Arguments ..
57 : INTEGER N,NNZ
58 : INTEGER IRET
59 : * ..
60 : * .. Array Arguments ..
61 : * .. Matrices stockées en Morse
62 : INTEGER ROWPTR(N+1)
63 : INTEGER COLIND(NNZ)
64 : REAL*8 VAL(NNZ)
65 : * .. Vecteurs
66 : C INTEGER TMPCOL(N)
67 : C REAL*8 TMPVAL(N)
68 : REAL*8 DINV(N)
69 : *
70 : * .. Variables locales
71 : * .. Parameters
72 : REAL*8 ZERO ,PETIT ,ONE
73 : PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
74 : * ..
75 : C Nombre de pivots petits
76 : INTEGER NBPIVP
77 : C Nombre de pivots inférieurs à 0
78 : INTEGER NBPIVI
79 : REAL*8 VALPIV,VALDIA
80 : C***
81 : IRET=0
82 : NBPIVP=0
83 : NBPIVI=0
84 : C On boucle sur les lignes
85 : DO 1 IN=1,N
86 : IDEB=ROWPTR(IN)
87 : IFIN=ROWPTR(IN+1)-1
88 : ICOL=IDEB
89 : IF (IDEB.LE.IFIN) THEN
90 : C On cherche le terme Aii
91 : 11 CONTINUE
92 : IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
93 : ICOL=ICOL+1
94 : GOTO 11
95 : ENDIF
96 : C On ne l'a pas trouvé
97 : IF (COLIND(ICOL).NE.IN) THEN
98 : WRITE(IOIMP,*) 'diag.',IN,'inexistante'
99 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
100 : WRITE(IOIMP,*) 'est impossible.'
101 : IRET=-1
102 : GOTO 9999
103 : ELSE
104 : C On l'a trouvé
105 : C Dii <- Dii+Aii
106 : C Dii <- 1/Dii
107 : VALDIA=VAL(ICOL)
108 : VALPIV=DINV(IN)+VALDIA
109 : IF (VALPIV.EQ.ZERO) THEN
110 : WRITE(IOIMP,*) 'Pivot',IN,'nul'
111 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
112 : WRITE(IOIMP,*) 'est impossible.'
113 : IRET=-2
114 : GOTO 9999
115 : ELSE
116 : IF (ABS(VALPIV).LT.PETIT) THEN
117 : NBPIVP=NBPIVP+1
118 : ENDIF
119 : IF (VALPIV.LT.ZERO) THEN
120 : NBPIVI=NBPIVI+1
121 : ENDIF
122 : ENDIF
123 : VALPIV=ONE/VALPIV
124 : DINV(IN)=VALPIV
125 : JCOL=ICOL+1
126 : JFIN=IFIN
127 : IF (JCOL.LE.JFIN) THEN
128 : C On parcourt les valeurs Aij (j>i).
129 : C Pour chaque valeur, on regarde ligne j s'il y a une
130 : C valeur Aji correspondante.
131 : C Si oui, Djj <- Djj - AjiDiiAij
132 : C Sinon, rien !
133 : 12 CONTINUE
134 : JN=COLIND(JCOL)
135 : C Recherche de Aji
136 : KDEB=ROWPTR(JN)
137 : KFIN=ROWPTR(JN+1)-1
138 : KCOL=KDEB
139 : 121 CONTINUE
140 : IF (COLIND(KCOL).LT.IN.AND.KCOL.LT.KFIN) THEN
141 : KCOL=KCOL+1
142 : GOTO 121
143 : ENDIF
144 : C On a trouve Aji !
145 : IF (COLIND(KCOL).EQ.IN) THEN
146 : DINV(JN)=DINV(JN)-
147 : $ (VAL(KCOL)*VALPIV*VAL(JCOL))
148 : ENDIF
149 : C On boucle...
150 : JCOL=JCOL+1
151 : IF (JCOL.LE.JFIN) GOTO 12
152 : ENDIF
153 : ENDIF
154 : ELSE
155 : WRITE(IOIMP,*) 'Ligne',IN,'vide'
156 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
157 : WRITE(IOIMP,*) 'est impossible.'
158 : IRET=-3
159 : GOTO 9999
160 : ENDIF
161 : 1 CONTINUE
162 : *
163 : * Warning(s)
164 : *
165 : IF (NBPIVP.GT.0) THEN
166 : WRITE(IOIMP,*) 'WARNING !'
167 : WRITE(IOIMP,*) 'mdilun :',NBPIVP,' |pivot(s)|<',PETIT
168 : ENDIF
169 : *! IF (NBPIVI.GT.0) THEN
170 : *! WRITE(IOIMP,*) 'WARNING !'
171 : *! WRITE(IOIMP,*) 'mdilun :',NBPIVI,' pivot(s) < 0'
172 : *! ENDIF
173 : *
174 : * Normal termination
175 : *
176 : RETURN
177 : *
178 : * Error handling
179 : *
180 : 9999 CONTINUE
181 : WRITE(IOIMP,*) 'An error was detected in mdilun.eso'
182 : RETURN
183 : *
184 : * End of MDILUN
185 : *
186 : END
187 :
188 :
189 :
190 :
191 :
192 :
193 :
194 :
195 :
© Cast3M 2003 - All rights reserved.
Disclaimer