1 : C PRMILU SOURCE MAGN 09/07/30 21:20:42 6458
2 : SUBROUTINE PRMILU(KMORS,KISA,MATRIK,RXMILU,IMPR,IRET)
3 : IMPLICIT INTEGER(I-N)
4 : IMPLICIT REAL*8 (A-H,O-Z)
5 : C***********************************************************************
6 : C NOM : PRMILU
7 : C DESCRIPTION :
8 : C Calcul du préconditionneur MILU(0) d'une matrice Morse.
9 : C MILU(0) : Modified Incomplete LU factorization of level 0
10 : C appelée aussi Choleski ou Crout incomplet modifié !
11 : C Ici, c'est meme un MILU(0) relaxé (avec RXMILU)...
12 : C
13 : C Le préconditionneur est une matrice stockée
14 : C au format MSR (Modified Sparse Row, stockage de l'inverse de la
15 : C diagonale) de meme profil que la matrice Morse (format CSR) qu'il
16 : C préconditionne.
17 : C Le profil et les valeurs du préconditionneur sont
18 : C stockés dans KIDMAT(6 et 7) (réutilisation de l'existant).
19 : C
20 : C Ce sous-programme est en fait une interface à :
21 : C memilu
22 : C qui est en Fortran presque pur (pour raison de rapidité)
23 : C et effectue la construction proprement dite du
24 : C préconditionneur.
25 : C
26 : C ATTENTION : pour une matrice A quelconque, la factorisation
27 : C --------- MILU(0) peut ne pas exister (pivot nul) ou avoir
28 : C des pivots négatifs MEME SI la factorisation
29 : C complète de A existe et n'a que des pivots
30 : C positifs.
31 : C
32 : C
33 : C LANGAGE : ESOPE
34 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
35 : C mél : gounand@semt2.smts.cea.fr
36 : C REFERENCE (bibtex-like) :
37 : C @BOOK{templates,
38 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
39 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
40 : C H. Van der Vorst},
41 : C TITLE={Templates for the Solution of Linear Systems :
42 : C Building Blocks for Iterative Methods},
43 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
44 : C -> URL : http://www.netlib.org/templates/Templates.html
45 : C***********************************************************************
46 : C APPELES : MEILU0
47 : C APPELES (E/S) : ECMORS
48 : C APPELE PAR : KRES2
49 : C***********************************************************************
50 : C ENTREES : MATRIK, RXMILU, IMPR
51 : C ENTREES/SORTIES : -
52 : C SORTIES : ILUM, ILUI (KIDMAT(6-7) dans MATRIK), IRET
53 : C CODE RETOUR (IRET) : 0 si ok
54 : C <0 si problème
55 : C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
56 : C on pioche dedans les informations nécessaires
57 : C (différents pointeurs, nb. de ddl...)
58 : C RXMILU : type REAL*8.
59 : C facteur de relaxation (normalement compris
60 : C entre 0.D0 et 1.D0)
61 : C 0.D0 : on se ramène à Choleski incomplet
62 : C 1.D0 : on se ramène à MILU(0) non relaxé.
63 : C IMPR : niveau d'impression
64 : C ILUM : pointeur sur segment PMORS de l'include SMMATRIK
65 : C profil morse du préconditionneur ILU(0)
66 : C =KIDMAT(6)=KMORS=KIDMAT(4) dans MATRIK
67 : C ILUI : pointeur sur segment IZA de l'include SMMATRIK
68 : C valeur du préconditionneur ILU(0)
69 : C =KIDMAT(7) dans MATRIK
70 : C***********************************************************************
71 : C VERSION : v1, 01/04/98, version initiale
72 : C HISTORIQUE : v1, 01/04/98, création
73 : C HISTORIQUE : 09/02/98, on ne construit pas le préconditionneur s'il
74 : C existe déjà.
75 : C HISTORIQUE : 20/12/99, interfaçage avec le nouveau memilu
76 : C Le préconditionneur est stocké au format MSR (Modified Sparse Row)
77 : C (voir la doc de Sparskit version 2+ (Youcef Saad))
78 : C HISTORIQUE :
79 : C HISTORIQUE :
80 : C***********************************************************************
81 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
82 : C en cas de modification de ce sous-programme afin de faciliter
83 : C la maintenance !
84 : C***********************************************************************
85 : -INC CCOPTIO
86 : -INC SMMATRIK
87 : -INC SMLENTI
88 : POINTEUR KMORS.PMORS
89 : POINTEUR KISA.IZA
90 : POINTEUR ILUM.PMORS
91 : POINTEUR ILUI.IZA
92 : POINTEUR IWORK.MLENTI
93 : REAL*8 RXMILU
94 : C***
95 : IRET=0
96 : IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans prmilu'
97 : C On récupère les segments utiles
98 : SEGACT MATRIK
99 : ILUM=KIDMAT(6)
100 : ILUI=KIDMAT(7)
101 : SEGDES MATRIK
102 : C Le préconditionneur est-il déjà construit ?
103 : IF ((ILUM.EQ.0).OR.(ILUI.EQ.0)) THEN
104 : C La matrice Morse et son préconditionneur ont le meme
105 : C profil
106 : SEGACT KMORS
107 : N=KMORS.IA(/1)-1
108 : NNZ=KMORS.JA(/1)
109 : SEGACT KISA
110 : NTT=N-1
111 : NJA=NNZ+1
112 : SEGINI ILUM
113 : NBVA=NNZ+1
114 : SEGINI ILUI
115 : JG=N
116 : SEGINI IWORK
117 : C Les boucles sont en Fortran pur
117 : CALL MEMILU(N,NNZ,KISA.A,KMORS.JA,KMORS.IA,
118 : $ ILUI.A,ILUM.JA,ILUM.IA,
119 : $ IWORK.LECT,
120 : $ RXMILU,
121 : $ IMPR,IRET)
122 : SEGSUP IWORK
123 : SEGDES ILUI
124 : SEGDES ILUM
125 : SEGDES KISA
126 : SEGDES KMORS
127 : IF (IRET.NE.0) GOTO 9999
128 : C
129 : C On stocke la factorisation obtenue du préconditionneur
130 : C
131 : SEGACT MATRIK*MOD
132 : KIDMAT(6)=ILUM
133 : KIDMAT(7)=ILUI
134 : SEGDES MATRIK
135 : C
136 : IF (IMPR.GT.6) THEN
137 : WRITE(IOIMP,*) 'création de la matrice Morse'
138 : WRITE(IOIMP,*) 'de pointeurs',ILUM,'et',ILUI
139 : IF (IMPR.GT.8) THEN
139 : CALL ECMORS(ILUM,ILUI,(IMPR-1))
140 : ENDIF
141 : ENDIF
142 : ELSE
143 : IF (IMPR.GT.6) THEN
144 : WRITE(IOIMP,*) 'Le préconditionneur est déjà construit :'
145 : WRITE(IOIMP,*) 'ILUM=',ILUM,' et ILUI=',ILUI
146 : IF (IMPR.GT.8) THEN
146 : CALL ECMORS(ILUM,ILUI,(IMPR-1))
147 : ENDIF
148 : ENDIF
149 : ENDIF
150 : *
151 : * Normal termination
152 : *
153 : RETURN
154 : *
155 : * Format handling
156 : *
157 : *
158 : * Error handling
159 : *
160 : 9999 CONTINUE
161 : WRITE(IOIMP,*) 'An error was detected in prmilu.eso'
162 : RETURN
163 : *
164 : * End of PRMILU
165 : *
166 : END
167 :
168 :
169 :
170 :
171 :
172 :
173 :
© Cast3M 2003 - All rights reserved.
Disclaimer