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