1 : C CGDILU SOURCE CHAT 05/01/12 21:53:22 5004
2 : SUBROUTINE CGDILU(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ DINV,
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 : CGDILU
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 la factorisation D-ILU de A.
16 : C
17 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
18 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
19 : C mél : gounand@semt2.smts.cea.fr
20 : C REFERENCE (bibtex-like) :
21 : C @BOOK{templates,
22 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
23 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
24 : C H. Van der Vorst},
25 : C TITLE={Templates for the Solution of Linear Systems :
26 : C Building Blocks for Iterative Methods},
27 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
28 : C -> URL : http://www.netlib.org/templates/Templates.html
29 : C***********************************************************************
30 : C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines)
31 : C DNRM2, DDOT (functions)
32 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
33 : C PSDILU : Preconditionner solve
34 : C (Montée-Descente adaptée à D-ILU)
35 : C***********************************************************************
36 : C ENTREES : N, NNZ,
37 : C ROWPTR, COLIND, VAL,
38 : C B, DINV, IMPR
39 : C ENTREES/SORTIES : X, R, P, Q, Z, ITER, RESID
40 : C SORTIES : LRES, IRET
41 : C CODE RETOUR (IRET) : 0 si ok
42 : C <0 si problème (non utilisé !)
43 : C >0 si l'algorithme n'a pas convergé
44 : C N : nombre de degrés de liberté
45 : C NNZ : nombre de valeurs non nulles de la matrice Morse
46 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
47 : C et valeurs de la matrice Morse
48 : C B : vecteur second membre
49 : C DINV : vecteur contenant les inverses des pivots
50 : C de la factorisation D-ILU de A.
51 : C IMPR : niveau d'impression (0..4)
52 : C X : vecteur des inconnues.
53 : C En entrée, contient l'estimation x(0).
54 : C En sortie, la solution trouvée (si la méthode
55 : C a convergé).
56 : C R,P,Q,Z: vecteurs de travail initialisés et détruits
57 : C dans prcg (on peut en utiliser moins ?).
58 : C ITER : type INTEGER.
59 : C En entrée, contient le nombre maximum
60 : C d'itérations à effectuer.
61 : C En sortie, contient le nombre d'itérations
62 : C réellement effectuées.
63 : C RESID : type REAL*8.
64 : C En entrée, la valeur maximale du résidu à
65 : C convergence de l'algorithme mesurée comme suit :
66 : C norme[L2](b - A*x) / norme[L2]( b )
67 : C En sortie, la valeur effective de cette mesure.
68 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
69 : C contient l'historique de la valeur de RESID en
70 : C fonction du nombre d'itérations effectuées.
71 : C
72 : C***********************************************************************
73 : C VERSION : v1, 02/04/98, version initiale
74 : C HISTORIQUE : v1, 02/04/98, création
75 : C HISTORIQUE :
76 : C HISTORIQUE :
77 : C***********************************************************************
78 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
79 : C en cas de modification de ce sous-programme afin de faciliter
80 : C la maintenance !
81 : C***********************************************************************
82 : *
83 : * .. Includes .. a commenter si mise au point ok
84 : -INC CCOPTIO
85 : INTEGER IMPR
86 : * .. Scalar Arguments ..
87 : INTEGER N,NNZ
88 : INTEGER ITER,IRET
89 : REAL*8 RESID
90 : * ..
91 : * .. Array Arguments ..
92 : * .. Matrices stockées en Morse
93 : INTEGER ROWPTR(N+1)
94 : INTEGER COLIND(NNZ)
95 : REAL*8 VAL(NNZ)
96 : * .. Vecteurs
97 : REAL*8 B(N)
98 : REAL*8 X(N)
99 : REAL*8 DINV(N)
100 : REAL*8 LRES(ITER)
101 : REAL*8 R(N),P(N),Q(N),Z(N)
102 : *
103 : *
104 : * .. Variables locales
105 : * .. Parameters
106 : REAL*8 ZERO , ONE
107 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
108 : * ..
109 : INTEGER MAXIT
110 : REAL*8 ALPHA,BETA,RHO,RHO1,TOL,BNRM2
111 : REAL*8 DDOT,DNRM2
112 : * ..
113 : * .. External subroutines and functions..
114 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2
115 : * EXTERNAL STOPTEST2
116 : * ..
117 : * .. Executable Statements ..
118 : *
119 :
120 : IRET = 0
121 : MAXIT = ITER
122 : TOL = RESID
123 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgdilu.eso'
124 : *
125 : * Calcul du résidu initial
126 : *
127 : C r(0)=b
127 : CALL DCOPY(N,B,1,R,1)
128 : IF(IMPR.GT.6) THEN
129 : WRITE(IOIMP,*) 'r(0)=b'
130 : ENDIF
131 : IF (DNRM2(N,X,1).NE.ZERO) THEN
132 : C r(0)=b-Ax
132 : CALL DMOMV ('No transpose',N,-ONE,
133 : $ NNZ,ROWPTR,COLIND,VAL,
134 : $ X,ONE,R)
135 : IF(IMPR.GT.6) THEN
136 : WRITE(IOIMP,*) 'r(0)=b-Ax'
137 : ENDIF
138 : * Je ne comprends pas cette ligne de code qui
139 : * vient de CGREVCOM.f
140 : * IF (DNRM2(N,R,1).LT.TOL) GOTO 30
141 : ENDIF
142 : *
143 : ITER = 0
144 :
145 : * On ne peut pas normer par B à cause des 1.E30
146 : BNRM2 = DNRM2(N,B,1)
147 : IF(IMPR.GT.5) THEN
148 : WRITE(IOIMP,*) '||B||=',BNRM2
149 : ENDIF
150 :
151 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
152 : RESID = DNRM2(N,R,1) / BNRM2
153 : * RESID = DNRM2(N,R,1)
154 : IF(IMPR.GT.5) THEN
155 : WRITE(IOIMP,*) 'Résidu initial=',RESID
156 : ENDIF
157 :
158 : IF (RESID.LE.TOL) THEN
159 : IRET = 0
160 : GOTO 30
161 : ENDIF
162 :
163 : *
164 : 10 CONTINUE
165 : *
166 : * Perform (Preconditioned) Conjugate Gradient Iteration.
167 : *
168 : ITER = ITER +1
169 : IF(IMPR.GT.6) THEN
170 : WRITE(IOIMP,*) 'ITER=',ITER
171 : ENDIF
172 :
173 : *
174 : * Preconditioner Solve (D-ILU)
175 : *
176 : C Mz(i-1)=r(i-1)
176 : CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
177 : $ DINV,
178 : $ Z,R)
179 : *
180 : *
181 : C rho(i-1)=r(i-1)'z(i-1)
182 : RHO = DDOT(N,R,1,Z,1)
183 : *
184 : * Compute direction vector P
185 : *
186 : IF (ITER.GT.1) THEN
187 : C beta(i-1)=rho(i-1)/rho(i-2)
188 : BETA = RHO / RHO1
189 : C p(i)=z(i-1)+beta(i-1)p(i-1)
189 : CALL DAXPY(N,BETA,P,1,Z,1)
189 : CALL DCOPY(N,Z,1,P,1)
190 : ELSE
191 : C p(1)=z(0)
191 : CALL DCOPY(N,Z,1,P,1)
192 : ENDIF
193 : *
194 : * Compute scalar ALPHA (save A*p to q)
195 : *
196 : C q(i)=Ap(i)
196 : CALL DMOMV('No transpose',N,ONE,
197 : $ NNZ,ROWPTR,COLIND,VAL,
198 : $ P,ZERO,Q)
199 : C alpha(i)=rho(i-1)/(p(i)'q(i))
200 : ALPHA = RHO / DDOT(N,P,1,Q,1)
201 : *
202 : * Compute current solution vector x
203 : *
204 : C x(i)=x(i-1)+alpha(i)p(i)
204 : CALL DAXPY(N,ALPHA,P,1,X,1)
205 : IF(IMPR.GT.7) THEN
206 : WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
207 : ENDIF
208 :
209 : *
210 : * Compute residual vector R, find norm,
211 : * then check for tolerance
212 : *
213 : C r(i)=r(i-1)-alpha(i)r(i)
213 : CALL DAXPY(N,-ALPHA,Q,1,R,1)
214 : *
215 : RESID = DNRM2(N,R,1) / BNRM2
216 : * RESID = DNRM2(N,R,1)
217 : IF(IMPR.GT.6) THEN
218 : WRITE(IOIMP,*) 'Résidu=',RESID
219 : ENDIF
220 : LRES(ITER)=RESID
221 : *
222 : * Iteration successful
223 : *
224 : IF (RESID.LE.TOL) THEN
225 : IRET = 0
226 : *
227 : * Iteration fails
228 : *
229 : ELSE IF (ITER.EQ.MAXIT) THEN
230 : IRET = 1
231 : *
232 : * Do some more
233 : *
234 : ELSE
235 : RHO1 = RHO
236 : GOTO 10
237 : ENDIF
238 : *
239 : * Normal termination
240 : *
241 : 30 CONTINUE
242 : RETURN
243 : *
244 : * End of CGDILU.
245 : *
246 : END
247 :
248 :
249 :
250 :
© Cast3M 2003 - All rights reserved.
Disclaimer