1 : C CGD SOURCE CHAT 05/01/12 21:53:19 5004
2 : SUBROUTINE CGD(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ MINV,
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 : CGD
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 diagonale de A (Jacobi).
16 : C
17 : C LANGAGE : FORTRAN 77 + chouilla 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***********************************************************************
34 : C ENTREES : N, NNZ,
35 : C ROWPTR, COLIND, VAL,
36 : C B, MINV, IMPR
37 : C ENTREES/SORTIES : X, R, P, Q, Z, ITER, RESID
38 : C SORTIES : LRES, IRET
39 : C CODE RETOUR (IRET) : 0 si ok
40 : C <0 si problème (non utilisé !)
41 : C >0 si l'algorithme n'a pas convergé
42 : C N : nombre de degrés de liberté
43 : C NNZ : nombre de valeurs non nulles de la matrice Morse
44 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
45 : C et valeurs de la matrice Morse
46 : C B : vecteur second membre
47 : C MINV : vecteur contenant les inverses des termes
48 : C diagonaux de A
49 : C IMPR : niveau d'impression (0..4)
50 : C X : vecteur des inconnues.
51 : C En entrée, contient l'estimation x(0).
52 : C En sortie, la solution trouvée (si la méthode
53 : C a convergé).
54 : C R,P,Q,Z: vecteurs de travail initialisés et détruits
55 : C dans prcg (on peut en utiliser moins ?).
56 : C ITER : type INTEGER.
57 : C En entrée, contient le nombre maximum
58 : C d'itérations à effectuer.
59 : C En sortie, contient le nombre d'itérations
60 : C réellement effectuées.
61 : C RESID : type REAL*8.
62 : C En entrée, la valeur maximale du résidu à
63 : C convergence de l'algorithme mesurée comme suit :
64 : C norme[L2](b - A*x) / norme[L2]( b )
65 : C En sortie, la valeur effective de cette mesure.
66 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
67 : C contient l'historique de la valeur de RESID en
68 : C fonction du nombre d'itérations effectuées.
69 : C
70 : C***********************************************************************
71 : C VERSION : v1, 02/04/98, version initiale
72 : C HISTORIQUE : v1, 02/04/98, création
73 : C HISTORIQUE :
74 : C HISTORIQUE :
75 : C***********************************************************************
76 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
77 : C en cas de modification de ce sous-programme afin de faciliter
78 : C la maintenance !
79 : C***********************************************************************
80 : *
81 : * .. Includes .. a commenter si mise au point ok
82 : -INC CCOPTIO
83 : INTEGER IMPR
84 : * .. Scalar Arguments ..
85 : INTEGER N,NNZ
86 : INTEGER ITER,IRET
87 : REAL*8 RESID
88 : * ..
89 : * .. Array Arguments ..
90 : * .. Matrices stockées en Morse
91 : INTEGER ROWPTR(N+1)
92 : INTEGER COLIND(NNZ)
93 : REAL*8 VAL(NNZ)
94 : * .. Vecteurs
95 : REAL*8 B(N)
96 : REAL*8 X(N)
97 : REAL*8 MINV(N)
98 : REAL*8 LRES(ITER)
99 : REAL*8 R(N),P(N),Q(N),Z(N)
100 : *
101 : *
102 : * .. Variables locales
103 : * .. Parameters
104 : REAL*8 ZERO , ONE
105 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
106 : * ..
107 : INTEGER MAXIT
108 : REAL*8 ALPHA,BETA,RHO,RHO1,TOL,BNRM2
109 : REAL*8 DDOT,DNRM2
110 : * ..
111 : * .. External subroutines and functions..
112 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2
113 : * EXTERNAL STOPTEST2
114 : * ..
115 : * .. Executable Statements ..
116 : *
117 :
118 : IRET = 0
119 : MAXIT = ITER
120 : TOL = RESID
121 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgd.eso'
122 : *
123 : * Calcul du résidu initial
124 : *
125 : C r(0)=b
125 : CALL DCOPY(N,B,1,R,1)
126 : IF(IMPR.GT.6) THEN
127 : WRITE(IOIMP,*) 'r(0)=b'
128 : ENDIF
129 : IF (DNRM2(N,X,1).NE.ZERO) THEN
130 : C r(0)=b-Ax
130 : CALL DMOMV ('No transpose',N,-ONE,
131 : $ NNZ,ROWPTR,COLIND,VAL,
132 : $ X,ONE,R)
133 : IF(IMPR.GT.6) THEN
134 : WRITE(IOIMP,*) 'r(0)=b-Ax'
135 : ENDIF
136 : * Je ne comprends pas cette ligne de code qui
137 : * vient de CGREVCOM.f
138 : * IF (DNRM2(N,R,1).LT.TOL) GOTO 30
139 : ENDIF
140 : *
141 : ITER = 0
142 :
143 : * On ne peut pas normer par B à cause des 1.E30
144 : BNRM2 = DNRM2(N,B,1)
145 : IF(IMPR.GT.5) THEN
146 : WRITE(IOIMP,*) '||B||=',BNRM2
147 : ENDIF
148 :
149 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
150 : RESID = DNRM2(N,R,1) / BNRM2
151 : * RESID = DNRM2(N,R,1)
152 : IF(IMPR.GT.5) THEN
153 : WRITE(IOIMP,*) 'Résidu initial=',RESID
154 : ENDIF
155 :
156 : IF (RESID.LE.TOL) THEN
157 : IRET = 0
158 : GOTO 30
159 : ENDIF
160 :
161 : *
162 : 10 CONTINUE
163 : *
164 : * Perform (Preconditioned) Conjugate Gradient Iteration.
165 : *
166 : ITER = ITER +1
167 : IF(IMPR.GT.6) THEN
168 : WRITE(IOIMP,*) 'ITER=',ITER
169 : ENDIF
170 :
171 : *
172 : * Preconditioner Solve
173 : *
174 : C Mz(i-1)=r(i-1)
175 : DO 100 IN=1,N
176 : Z(IN)=R(IN)*MINV(IN)
177 : 100 CONTINUE
178 : *
179 : *
180 : C rho(i-1)=r(i-1)'z(i-1)
181 : RHO = DDOT(N,R,1,Z,1)
182 : *
183 : * Compute direction vector P
184 : *
185 : IF (ITER.GT.1) THEN
186 : C beta(i-1)=rho(i-1)/rho(i-2)
187 : BETA = RHO / RHO1
188 : C p(i)=z(i-1)+beta(i-1)p(i-1)
188 : CALL DAXPY(N,BETA,P,1,Z,1)
188 : CALL DCOPY(N,Z,1,P,1)
189 : ELSE
190 : C p(1)=z(0)
190 : CALL DCOPY(N,Z,1,P,1)
191 : ENDIF
192 : *
193 : * Compute scalar ALPHA (save A*p to q)
194 : *
195 : C q(i)=Ap(i)
195 : CALL DMOMV ('No transpose',N,ONE,
196 : $ NNZ,ROWPTR,COLIND,VAL,
197 : $ P,ZERO,Q)
198 : C alpha(i)=rho(i-1)/(p(i)'q(i))
199 : ALPHA = RHO / DDOT(N,P,1,Q,1)
200 : *
201 : * Compute current solution vector x
202 : *
203 : C x(i)=x(i-1)+alpha(i)p(i)
203 : CALL DAXPY(N,ALPHA,P,1,X,1)
204 : IF(IMPR.GT.7) THEN
205 : WRITE(IOIMP,*) 'X=',(X(III),III=1,N)
206 : ENDIF
207 :
208 : *
209 : * Compute residual vector R, find norm,
210 : * then check for tolerance
211 : *
212 : C r(i)=r(i-1)-alpha(i)r(i)
212 : CALL DAXPY(N,-ALPHA,Q,1,R,1)
213 : *
214 : RESID = DNRM2(N,R,1) / BNRM2
215 : * RESID = DNRM2(N,R,1)
216 : IF(IMPR.GT.6) THEN
217 : WRITE(IOIMP,*) 'Résidu=',RESID
218 : ENDIF
219 : LRES(ITER)=RESID
220 :
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 CGD.
245 : *
246 : END
247 :
248 :
249 :
250 :
© Cast3M 2003 - All rights reserved.
Disclaimer