1 : C CGILU0 SOURCE CHAT 05/01/12 21:53:28 5004
2 : SUBROUTINE CGILU0(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ JU,JLU,ALU,
5 : $ LRES,
6 : $ R,P,Q,Z,
7 : $ ITER,RESID,IMPR,IRET)
8 : IMPLICIT INTEGER(I-N)
9 : IMPLICIT REAL*8 (A-H,O-Z)
10 : C***********************************************************************
11 : C NOM : CGILU0
12 : C DESCRIPTION :
13 : C Résolution d'un système linéaire Ax=b
14 : C par la méthode du gradient conjugué
15 : C préconditionné par une factorisation ILU(0) (Crout incomplet)
16 : C ou MILU(0) (Crout inc. modifié relaxé) de A.
17 : C
18 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les 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 (BLAS 1) : DCOPY, DAXPY (subroutines)
32 : C DNRM2, DDOT (functions)
33 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
34 : C PSILU0 : Preconditionner solve
35 : C (Montée-Descente adaptée à ILU(0))
36 : C***********************************************************************
37 : C ENTREES : N, NNZ,
38 : C ROWPTR, COLIND, VAL,
39 : C B, JU, JLU, ALU, IMPR
40 : C ENTREES/SORTIES : X, R, P, Q, Z, ITER, RESID
41 : C SORTIES : LRES, IRET
42 : C CODE RETOUR (IRET) : 0 si ok
43 : C <0 si problème (non utilisé !)
44 : C >0 si l'algorithme n'a pas convergé
45 : C N : nombre de degrés de liberté
46 : C NNZ : nombre de valeurs non nulles de la matrice Morse
47 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
48 : C et valeurs de la matrice Morse A
49 : C B : vecteur second membre
50 : C ILUVAL : valeurs de la matrice de la factorisation
51 : C ILU(0) ou MILU(0) de A.
52 : C Cette matrice est stockée en Morse, elle a le
53 : C meme profil que A.
54 : C IMPR : niveau d'impression (0..4)
55 : C X : vecteur des inconnues.
56 : C En entrée, contient l'estimation x(0).
57 : C En sortie, la solution trouvée (si la méthode
58 : C a convergé).
59 : C R,P,Q,Z: vecteurs de travail initialisés et détruits
60 : C dans prcg (on peut en utiliser moins ?).
61 : C ITER : type INTEGER.
62 : C En entrée, contient le nombre maximum
63 : C d'itérations à effectuer.
64 : C En sortie, contient le nombre d'itérations
65 : C réellement effectuées.
66 : C RESID : type REAL*8.
67 : C En entrée, la valeur maximale du résidu à
68 : C convergence de l'algorithme mesurée comme suit :
69 : C norme[L2](b - A*x) / norme[L2]( b )
70 : C En sortie, la valeur effective de cette mesure.
71 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
72 : C contient l'historique de la valeur de RESID en
73 : C fonction du nombre d'itérations effectuées.
74 : C
75 : C***********************************************************************
76 : C VERSION : 20/12/99
77 : C HISTORIQUE : v1, 02/04/98, création
78 : C HISTORIQUE : 20/12/99
79 : C Interfaçage avec le nouveau psilu0 (factorisations incomplètes
80 : C stockées au format MSR (Modified Sparse Row).
81 : C HISTORIQUE :
82 : C HISTORIQUE :
83 : C***********************************************************************
84 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
85 : C en cas de modification de ce sous-programme afin de faciliter
86 : C la maintenance !
87 : C***********************************************************************
88 : *
89 : * .. Includes .. a commenter si mise au point ok
90 : -INC CCOPTIO
91 : INTEGER IMPR
92 : * .. Scalar Arguments ..
93 : INTEGER N,NNZ
94 : INTEGER ITER,IRET
95 : REAL*8 RESID
96 : * ..
97 : * .. Array Arguments ..
98 : * .. Matrices stockées en Morse
99 : INTEGER ROWPTR(N+1)
100 : INTEGER COLIND(NNZ)
101 : REAL*8 VAL(NNZ)
102 : INTEGER JU(N)
103 : INTEGER JLU(NNZ+1)
104 : REAL*8 ALU(NNZ)
105 : * .. Vecteurs
106 : REAL*8 B(N)
107 : REAL*8 X(N)
108 : REAL*8 LRES(ITER)
109 : REAL*8 R(N),P(N),Q(N),Z(N)
110 : *
111 : *
112 : * .. Variables locales
113 : * .. Parameters
114 : REAL*8 ZERO , ONE
115 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
116 : * ..
117 : INTEGER MAXIT,III
118 : REAL*8 ALPHA,BETA,RHO,RHO1,TOL,BNRM2
119 : REAL*8 DDOT,DNRM2
120 : * ..
121 : * .. External subroutines and functions..
122 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2
123 : * EXTERNAL STOPTEST2
124 : * ..
125 : * .. Executable Statements ..
126 : *
127 :
128 : IRET = 0
129 : MAXIT = ITER
130 : TOL = RESID
131 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgilu0.eso'
132 : *
133 : * Calcul du résidu initial
134 : *
135 : C r(0)=b
135 : CALL DCOPY(N,B,1,R,1)
136 : IF(IMPR.GT.7) THEN
137 : WRITE(IOIMP,*) 'r(0)=b'
138 : ENDIF
139 : IF (DNRM2(N,X,1).NE.ZERO) THEN
140 : C r(0)=b-Ax
140 : CALL DMOMV ('No transpose',N,-ONE,
141 : $ NNZ,ROWPTR,COLIND,VAL,
142 : $ X,ONE,R)
143 : IF(IMPR.GT.7) THEN
144 : WRITE(IOIMP,*) 'r(0)=b-Ax'
145 : ENDIF
146 : * Je ne comprends pas cette ligne de code qui
147 : * vient de CGREVCOM.f
148 : * IF (DNRM2(N,R,1).LT.TOL) GOTO 30
149 : ENDIF
150 : *
151 : ITER = 0
152 :
153 : * On ne peut pas normer par B à cause des 1.E30
154 : BNRM2 = DNRM2(N,B,1)
155 : IF(IMPR.GT.5) THEN
156 : WRITE(IOIMP,*) '||B||=',BNRM2
157 : ENDIF
158 :
159 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
160 : RESID = DNRM2(N,R,1) / BNRM2
161 : * RESID = DNRM2(N,R,1)
162 : IF(IMPR.GT.5) THEN
163 : WRITE(IOIMP,*) 'Résidu initial=',RESID
164 : ENDIF
165 :
166 : IF (RESID.LE.TOL) THEN
167 : IRET = 0
168 : GOTO 30
169 : ENDIF
170 :
171 : *
172 : 10 CONTINUE
173 : *
174 : * Perform (Preconditioned) Conjugate Gradient Iteration.
175 : *
176 : ITER = ITER +1
177 : IF(IMPR.GT.7) THEN
178 : WRITE(IOIMP,*) 'ITER=',ITER
179 : ENDIF
180 :
181 : *
182 : * Preconditioner Solve (ILU(0))
183 : *
184 : C Mz(i-1)=r(i-1)
184 : CALL PSILU0(N,NNZ,JU,JLU,ALU,
185 : $ Z,R)
186 : *
187 : C rho(i-1)=r(i-1)'z(i-1)
188 : RHO = DDOT(N,R,1,Z,1)
189 : *
190 : * Compute direction vector P
191 : *
192 : IF (ITER.GT.1) THEN
193 : C beta(i-1)=rho(i-1)/rho(i-2)
194 : BETA = RHO / RHO1
195 : C p(i)=z(i-1)+beta(i-1)p(i-1)
195 : CALL DAXPY(N,BETA,P,1,Z,1)
195 : CALL DCOPY(N,Z,1,P,1)
196 : ELSE
197 : C p(1)=z(0)
197 : CALL DCOPY(N,Z,1,P,1)
198 : ENDIF
199 : *
200 : * Compute scalar ALPHA (save A*p to q)
201 : *
202 : C q(i)=Ap(i)
202 : CALL DMOMV('No transpose',N,ONE,
203 : $ NNZ,ROWPTR,COLIND,VAL,
204 : $ P,ZERO,Q)
205 : C alpha(i)=rho(i-1)/(p(i)'q(i))
206 : ALPHA = RHO / DDOT(N,P,1,Q,1)
207 : *
208 : * Compute current solution vector x
209 : *
210 : C x(i)=x(i-1)+alpha(i)p(i)
210 : CALL DAXPY(N,ALPHA,P,1,X,1)
211 : IF(IMPR.GT.8) THEN
212 : WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
213 : ENDIF
214 :
215 : *
216 : * Compute residual vector R, find norm,
217 : * then check for tolerance
218 : *
219 : C r(i)=r(i-1)-alpha(i)r(i)
219 : CALL DAXPY(N,-ALPHA,Q,1,R,1)
220 : *
221 : RESID = DNRM2(N,R,1) / BNRM2
222 : * RESID = DNRM2(N,R,1)
223 : IF(IMPR.GT.7) THEN
224 : WRITE(IOIMP,*) 'Résidu=',RESID
225 : ENDIF
226 : LRES(ITER)=RESID
227 : *
228 : * Iteration successful
229 : *
230 : IF (RESID.LE.TOL) THEN
231 : IRET = 0
232 : *
233 : * Iteration fails
234 : *
235 : ELSE IF (ITER.EQ.MAXIT) THEN
236 : IRET = 1
237 : *
238 : * Do some more
239 : *
240 : ELSE
241 : RHO1 = RHO
242 : GOTO 10
243 : ENDIF
244 : *
245 : * Normal termination
246 : *
247 : 30 CONTINUE
248 : RETURN
249 : *
250 : * End of CGILU0.
251 : *
252 : END
253 :
254 :
255 :
256 :
257 :
258 :
© Cast3M 2003 - All rights reserved.
Disclaimer