1 : C GMRD SOURCE CHAT 05/01/13 00:18:38 5004
2 : SUBROUTINE GMRD(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ MINV,
5 : $ LRES,
6 : $ R,W,AV,
7 : $ Y,GCOS,GSIN,S,
8 : $ V,H,
9 : $ ITER,RESTRT,RESID,IMPR,IRET)
10 : IMPLICIT INTEGER(I-N)
11 : IMPLICIT REAL*8 (A-H,O-Z)
12 : C***********************************************************************
13 : C NOM : GMRD
14 : C DESCRIPTION :
15 : C Résolution d'un système linéaire Ax=b
16 : C par une méthode GMRES(m) (restarted Generalized Minimal
17 : C Residual) préconditionné par la diagonale de A (Jacobi).
18 : C
19 : C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
20 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
21 : C mél : gounand@semt2.smts.cea.fr
22 : C REFERENCE (bibtex-like) :
23 : C @BOOK{templates,
24 : C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
25 : C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
26 : C H. Van der Vorst},
27 : C TITLE={Templates for the Solution of Linear Systems :
28 : C Building Blocks for Iterative Methods},
29 : C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
30 : C -> URL : http://www.netlib.org/templates/Templates.html
31 : C***********************************************************************
32 : C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines)
33 : C DNRM2, DDOT, DSCAL (functions)
34 : C APPELES (BLAS 2) : DTRSV (subroutines)
35 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
36 : C***********************************************************************
37 : C ENTREES : N, NNZ,
38 : C ROWPTR, COLIND, VAL,
39 : C B, MINV, RESTRT, IMPR
40 : C ENTREES/SORTIES : X,
41 : C R,W,AV,
42 : C Y,GCOS,GSIN,S,
43 : C V,H,
44 : C ITER, RESID
45 : C SORTIES : LRES, IRET
46 : C CODE RETOUR (IRET) : 0 si ok
47 : C <0 si problème (cf. BRTOL)
48 : C >0 si l'algorithme n'a pas convergé
49 : C N : nombre de degrés de liberté
50 : C NNZ : nombre de valeurs non nulles de la matrice Morse
51 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
52 : C et valeurs de la matrice Morse
53 : C B : vecteur second membre
54 : C MINV : vecteur contenant les inverses des termes
55 : C diagonaux de A
56 : C RESTRT : INTEGER
57 : C paramètre de redémarrage (restart) (i.e. m) pour
58 : C la méthode GMRES(m). (cf. prgmr.eso)
59 : C IMPR : niveau d'impression (0..4)
60 : C X : vecteur des inconnues.
61 : C En entrée, contient l'estimation x(0).
62 : C En sortie, la solution trouvée (si la méthode
63 : C a convergé).
64 : C R,W,AV : vecteurs de travail initialisés et détruits
65 : C Y,GCOS,GSIN,S dans prgmr.
66 : C V,H : tableaux de travail...
67 : C ITER : type INTEGER.
68 : C En entrée, contient le nombre maximum
69 : C d'itérations à effectuer.
70 : C En sortie, contient le nombre d'itérations
71 : C réellement effectuées.
72 : C RESID : type REAL*8.
73 : C En entrée, la valeur maximale du résidu à
74 : C convergence de l'algorithme mesurée comme suit :
75 : C norme[L2](b - A*x) / norme[L2]( b )
76 : C En sortie, la valeur effective de cette mesure.
77 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
78 : C contient l'historique de la valeur de RESID en
79 : C fonction du nombre d'itérations effectuées.
80 : C
81 : C***********************************************************************
82 : C VERSION : v1, 16/06/98, version initiale
83 : C HISTORIQUE : v1, 16/06/98, création
84 : C HISTORIQUE :
85 : C HISTORIQUE :
86 : C***********************************************************************
87 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
88 : C en cas de modification de ce sous-programme afin de faciliter
89 : C la maintenance !
90 : C***********************************************************************
91 : *
92 : * .. Includes .. a commenter si mise au point ok
93 : -INC CCOPTIO
94 : INTEGER IMPR
95 : * .. Scalar Arguments ..
96 : INTEGER N,NNZ
97 : INTEGER ITER,RESTRT,IRET
98 : REAL*8 RESID
99 : * ..
100 : * .. Array Arguments ..
101 : * .. Matrices stockées en Morse
102 : INTEGER ROWPTR(N+1)
103 : INTEGER COLIND(NNZ)
104 : REAL*8 VAL(NNZ)
105 : * .. Vecteurs
106 : REAL*8 B(N)
107 : REAL*8 X(N)
108 : REAL*8 MINV(N)
109 : REAL*8 LRES(ITER)
110 : * .. Tableaux de travail
111 : REAL*8 R(N),W(N),AV(N)
112 : REAL*8 Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT)
113 : REAL*8 S(RESTRT+1)
114 : REAL*8 V(N,RESTRT+1)
115 : REAL*8 H(RESTRT+1,RESTRT)
116 : *
117 : *
118 : * .. Variables locales
119 : * .. Parameters
120 : REAL*8 ZERO,ONE
121 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
122 : * ..
123 : INTEGER I,MAXIT,J,K,IN
124 : REAL*8 BNRM2,TOL,DDOT,DNRM2
125 : REAL*8 H1,H2,RAY
126 : * ..
127 : * .. External subroutines and functions..
128 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DTRSV
129 : EXTERNAL DMOMV
130 : * ..
131 : * .. Executable Statements ..
132 : *
133 : *
134 : IRET = 0
135 : MAXIT = ITER
136 : TOL = RESID
137 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmrd.eso'
138 : *
139 : * Set initial residual.
140 : *
141 : C r=b
141 : CALL DCOPY(N,B,1,R,1)
142 : IF (DNRM2(N,X,1).NE.ZERO) THEN
143 : C r=b-Ax
143 : CALL DMOMV ('No transpose',N,-ONE,
144 : $ NNZ,ROWPTR,COLIND,VAL,
145 : $ X,ONE,R)
146 : ENDIF
147 : *
148 : ITER = 0
149 :
150 : BNRM2 = DNRM2(N,B,1)
151 : IF(IMPR.GT.6) THEN
152 : WRITE(IOIMP,*) '||B||=',BNRM2
153 : ENDIF
154 :
155 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
156 : RESID = DNRM2(N,R,1) / BNRM2
157 : * RESID = DNRM2(N,R,1)
158 : IF(IMPR.GT.6) THEN
159 : WRITE(IOIMP,*) 'Résidu initial=',RESID
160 : ENDIF
161 :
162 : IF (RESID.LE.TOL) GOTO 30
163 : *
164 : 10 CONTINUE
165 : *
166 : * Construct the first column of V, and initialize S to the
167 : * elementary vector E1 scaled by RNORM.
168 : *
169 : C Preconditionner solve : Mv(1)=r
169 : C CALL DCOPY(N,R,1,V(1,1),1)
170 : DO 100 IN=1,N
171 : V(IN,1)=R(IN)*MINV(IN)
172 : 100 CONTINUE
173 :
174 : C s=||v(1)|| e(1)
175 : C v(1)=v(1) / ||v(1)||
176 : S(1) = DNRM2(N,V(1,1),1)
176 : CALL DSCAL(N,ONE/S(1),V(1,1),1)
177 : *
178 : DO 50 I = 1, RESTRT
179 : ITER = ITER + 1
180 : C av=A v(i)
180 : CALL DMOMV ('No transpose',N,ONE,
181 : $ NNZ,ROWPTR,COLIND,VAL,
182 : $ V(1,I),ZERO,AV)
183 : C Preconditionner solve : Mw=av
183 : C CALL DCOPY(N,AV,1,W,1)
184 : DO 200 IN=1,N
185 : W(IN)=AV(IN)*MINV(IN)
186 : 200 CONTINUE
187 : *
188 : * Construct I-th column of H so that it is orthonormal to
189 : * the previous I-1 columns.
190 : *
190 : C CALL ORTHOH(I,N,RESTRT,WORK2(1,I+H-1),WORK(1,V),WORK(1,W))
191 : DO 60 K = 1, I
192 : H(K,I) = DDOT(N,W,1,V(1,K),1)
192 : CALL DAXPY(N,-H(K,I),V(1,K),1,W,1)
193 : 60 CONTINUE
194 : H(I+1,I)=DNRM2(N,W,1)
194 : CALL DCOPY(N,W,1,V(1,I+1),1)
194 : CALL DSCAL(N,ONE/H(I+1,I),V(1,I+1),1)
195 : *
196 : C IF (I.GT.1)
197 : *
198 : * Apply Givens rotations to the I-th column of H. This
199 : * effectively reduces the Hessenberg matrix to upper
200 : * triangular form during the RESTRT iterations.
201 : *
201 : C $ CALL APPGIV(I,RESTRT,WORK2(1,I+H-1),WORK2(1,GIV))
202 : DO 70 K = 1, I-1
203 : H1= (GCOS(K)*H(K,I)) + (GSIN(K)*H(K+1,I))
204 : H2=-(GSIN(K)*H(K,I)) + (GCOS(K)*H(K+1,I))
205 : H(K,I) =H1
206 : H(K+1,I)=H2
207 : 70 CONTINUE
208 : *
209 : * Approximate residual norm. Check tolerance. If okay, compute
210 : * final approximation vector X and quit.
211 : *
212 : RAY=SQRT((H(I,I)*H(I,I)) + (H(I+1,I)*H(I+1,I)))
213 : GCOS(I) = H(I,I) / RAY
214 : GSIN(I) = H(I+1,I) / RAY
215 : H(I,I)=RAY
216 : H(I+1,I)=ZERO
217 : S(I+1)=-GSIN(I)*S(I)
218 : S(I) = GCOS(I)*S(I)
219 : RESID=ABS(S(I+1)) / BNRM2
220 : LRES(ITER)=RESID
221 : IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN
222 : * Solve the system Hy=S
223 : * DO 501 J=I,1,-1
224 : * Y(J)=S(J) / H(J,J)
225 : * DO 502 K=J-1,1,-1
226 : * S(K)=S(K)-H(K,J)*Y(J)
227 : * 502 CONTINUE
228 : * 501 CONTINUE
228 : CALL DCOPY(I,S,1,Y,1)
228 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H,RESTRT+1,Y,1)
229 : * Updating solution
229 : C CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
230 : C $ WORK(1,S),WORK(1,V))
231 : DO 503 J=I,1,-1
231 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
232 : 503 CONTINUE
233 : IF (RESID.LE.TOL) THEN
234 : GOTO 30
235 : ENDIF
236 : IF (ITER.EQ.MAXIT) THEN
237 : GOTO 40
238 : ENDIF
239 : ENDIF
240 : 50 CONTINUE
241 : *
242 : * Compute current solution vector X.
243 : *
244 : * Solve the system Hy=S
245 : * DO 301 J=RESTRT,1,-1
246 : * Y(J)=S(J) / H(J,J)
247 : * DO 302 K=J-1,1,-1
248 : * S(K)=S(K)-H(K,J)*Y(J)
249 : * 302 CONTINUE
250 : * 301 CONTINUE
250 : CALL DCOPY(RESTRT,S,1,Y,1)
250 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',
251 : $ RESTRT,H,RESTRT+1,Y,1)
252 : * Updating solution
252 : C CALL UPDATE(RESTRT,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
253 : C $ WORK(1,S),WORK(1,V))
254 : DO 303 J=RESTRT,1,-1
254 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
255 : 303 CONTINUE
256 : *
257 : * Compute residual vector R, find norm,
258 : * then check for tolerance.
259 : *
259 : CALL DCOPY(N,B,1,R,1)
259 : CALL DMOMV ('No transpose',N,-ONE,
260 : $ NNZ,ROWPTR,COLIND,VAL,
261 : $ X,ONE,R)
262 : S(RESTRT+1) = DNRM2(N,R,1)
263 : RESID = S(RESTRT+1) / BNRM2
264 : LRES(ITER)=RESID
265 : IF (RESID.LE.TOL) GOTO 30
266 : *
267 : * Iteration fails.
268 : *
269 : IF (ITER.EQ.MAXIT) THEN
270 : GOTO 40
271 : ENDIF
272 : *
273 : * Do some more
274 : *
275 : GOTO 10
276 : *
277 : * Iteration successful
278 : *
279 : 30 CONTINUE
280 : IRET = 0
281 : RETURN
282 : *
283 : * Iteration fails.
284 : *
285 : 40 CONTINUE
286 : IRET = 1
287 : RETURN
288 : *
289 : * End of GMRD.
290 : *
291 : END
292 :
293 :
294 :
295 :
296 :
297 :
298 :
© Cast3M 2003 - All rights reserved.
Disclaimer