1 : C MDILUS SOURCE CHAT 05/01/13 01:40:17 5004
2 : SUBROUTINE MDILUS(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 : MDILUS
10 : C DESCRIPTION :
11 : C Calcul du préconditionneur D-ILU d'une matrice Morse
12 : C SYMETRIQUE.
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 + chouhia 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 : REAL*8 DINV(N)
67 : *
68 : * .. Variables locales
69 : * .. Parameters
70 : REAL*8 ZERO ,PETIT ,ONE
71 : PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
72 : * ..
73 : C Nombre de pivots petits
74 : INTEGER NBPIVP
75 : C Nombre de pivots inférieurs à 0
76 : INTEGER NBPIVI
77 : REAL*8 VALPIV,VALDIA
78 : C***
79 : IRET=0
80 : NBPIVP=0
81 : NBPIVI=0
82 : C On boucle sur les lignes
83 : DO 1 IN=1,N
84 : IDEB=ROWPTR(IN)
85 : IFIN=ROWPTR(IN+1)-1
86 : ICOL=IDEB
87 : IF (IDEB.LE.IFIN) THEN
88 : C On cherche le terme Aii
89 : 11 CONTINUE
90 : IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
91 : ICOL=ICOL+1
92 : GOTO 11
93 : ENDIF
94 : C On ne l'a pas trouvé
95 : IF (COLIND(ICOL).NE.IN) THEN
96 : WRITE(IOIMP,*) 'diag.',IN,'inexistante'
97 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
98 : WRITE(IOIMP,*) 'est impossible.'
99 : IRET=-1
100 : GOTO 9999
101 : ENDIF
102 : C On l'a trouvé
103 : C Dii <- Dii+Aii
104 : C Dii <- 1/Dii
105 : VALDIA=VAL(ICOL)
106 : VALPIV=DINV(IN)+VALDIA
107 : IF (VALPIV.EQ.ZERO) THEN
108 : WRITE(IOIMP,*) 'Pivot',IN,'nul'
109 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
110 : WRITE(IOIMP,*) 'est impossible.'
111 : IRET=-2
112 : GOTO 9999
113 : ELSE
114 : IF (ABS(VALPIV).LT.PETIT) THEN
115 : NBPIVP=NBPIVP+1
116 : ENDIF
117 : IF (VALPIV.LT.ZERO) THEN
118 : NBPIVI=NBPIVI+1
119 : ENDIF
120 : ENDIF
121 : VALPIV=ONE/VALPIV
122 : DINV(IN)=VALPIV
123 : C pour j>i : Djj <- Djj-AjiDiiAij
124 : JCOL=ICOL+1
125 : JFIN=IFIN
126 : IF (JCOL.LE.JFIN) THEN
127 : 12 CONTINUE
128 : JN=COLIND(JCOL)
129 : AIJ=VAL(JCOL)
130 : DINV(JN)=DINV(JN)-(VALPIV*AIJ*AIJ)
131 : JCOL=JCOL+1
132 : IF (JCOL.LE.JFIN) GOTO 12
133 : ENDIF
134 : ELSE
135 : WRITE(IOIMP,*) 'Ligne',IN,'vide'
136 : WRITE(IOIMP,*) 'le préconditionnement D-ILU'
137 : WRITE(IOIMP,*) 'est impossible.'
138 : IRET=-3
139 : GOTO 9999
140 : ENDIF
141 : 1 CONTINUE
142 : *
143 : * Warning(s)
144 : *
145 : IF (NBPIVP.GT.0) THEN
146 : WRITE(IOIMP,*) 'WARNING !'
147 : WRITE(IOIMP,*) 'mdilus :',NBPIVP,' |pivot(s)|<',PETIT
148 : ENDIF
149 : IF (NBPIVI.GT.0) THEN
150 : WRITE(IOIMP,*) 'WARNING !'
151 : WRITE(IOIMP,*) 'mdilus :',NBPIVI,' pivot(s) < 0'
152 : ENDIF
153 : *
154 : * Normal termination
155 : *
156 : RETURN
157 : *
158 : * Format handling
159 : *
160 : *
161 : * Error handling
162 : *
163 : 9999 CONTINUE
164 : WRITE(IOIMP,*) 'An error was detected in mdilus.eso'
165 : RETURN
166 : *
167 : * End of MDILUS
168 : *
169 : END
170 :
171 :
172 :
173 :
174 :
175 :
176 :
177 :
© Cast3M 2003 - All rights reserved.
Disclaimer