1 : C GMR SOURCE CHAT 05/01/13 00:18:45 5004
2 : SUBROUTINE GMR(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ LRES,
5 : $ R,W,AV,
6 : $ Y,GCOS,GSIN,S,
7 : $ V,H,
8 : $ ITER,RESTRT,RESID,IMPR,IRET)
9 : IMPLICIT INTEGER(I-N)
10 : IMPLICIT REAL*8 (A-H,O-Z)
11 : C***********************************************************************
12 : C NOM : GMR
13 : C DESCRIPTION :
14 : C Résolution d'un système linéaire Ax=b
15 : C par une méthode GMRES(m) (restarted Generalized Minimal
16 : C Residual) non préconditionnée.
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, DSCAL (functions)
33 : C APPELES (BLAS 2) : DTRSV (subroutines)
34 : C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
35 : C***********************************************************************
36 : C ENTREES : N, NNZ,
37 : C ROWPTR, COLIND, VAL,
38 : C B, RESTRT, IMPR
39 : C ENTREES/SORTIES : X,
40 : C R,W,AV,
41 : C Y,GCOS,GSIN,S,
42 : C V,H,
43 : C ITER, RESID
44 : C SORTIES : LRES, IRET
45 : C CODE RETOUR (IRET) : 0 si ok
46 : C <0 si problème (cf. BRTOL)
47 : C >0 si l'algorithme n'a pas convergé
48 : C N : nombre de degrés de liberté
49 : C NNZ : nombre de valeurs non nulles de la matrice Morse
50 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
51 : C et valeurs de la matrice Morse
52 : C B : vecteur second membre
53 : C RESTRT : INTEGER
54 : C paramètre de redémarrage (restart) (i.e m) pour
55 : C la méthode GMRES(m). (cf. prgmr.eso)
56 : C IMPR : niveau d'impression (0..4)
57 : C X : vecteur des inconnues.
58 : C En entrée, contient l'estimation x(0).
59 : C En sortie, la solution trouvée (si la méthode
60 : C a convergé).
61 : C R,W,AV : vecteurs de travail initialisés et détruits
62 : C Y,GCOS,GSIN,S dans prgmr.
63 : C V,H : tableaux de travail...
64 : C ITER : type INTEGER.
65 : C En entrée, contient le nombre maximum
66 : C d'itérations à effectuer.
67 : C En sortie, contient le nombre d'itérations
68 : C réellement effectuées.
69 : C RESID : type REAL*8.
70 : C En entrée, la valeur maximale du résidu à
71 : C convergence de l'algorithme mesurée comme suit :
72 : C norme[L2](b - A*x) / norme[L2]( b )
73 : C En sortie, la valeur effective de cette mesure.
74 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
75 : C contient l'historique de la valeur de RESID en
76 : C fonction du nombre d'itérations effectuées.
77 : C
78 : C***********************************************************************
79 : C VERSION : v1, 16/06/98, version initiale
80 : C HISTORIQUE : v1, 16/06/98, création
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,RESTRT,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 : * .. Vecteurs
103 : REAL*8 B(N)
104 : REAL*8 X(N)
105 : REAL*8 LRES(ITER)
106 : * .. Tableaux de travail
107 : REAL*8 R(N),W(N),AV(N)
108 : REAL*8 Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT)
109 : REAL*8 S(RESTRT+1)
110 : REAL*8 V(N,RESTRT+1)
111 : REAL*8 H(RESTRT+1,RESTRT)
112 : *
113 : *
114 : * .. Variables locales
115 : * .. Parameters
116 : REAL*8 ZERO,ONE
117 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
118 : * ..
119 : INTEGER I,MAXIT,J,K
120 : REAL*8 BNRM2,TOL,DDOT,DNRM2
121 : REAL*8 H1,H2,RAY
122 : * ..
123 : * .. External subroutines and functions..
124 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DTRSV,DMOMV
125 : * ..
126 : * .. Executable Statements ..
127 : *
128 : *
129 : IRET = 0
130 : MAXIT = ITER
131 : TOL = RESID
132 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmr.eso'
133 : *
134 : * Set initial residual.
135 : *
136 : C r=b
136 : CALL DCOPY(N,B,1,R,1)
137 : IF (DNRM2(N,X,1).NE.ZERO) THEN
138 : C r=b-Ax
138 : CALL DMOMV ('No transpose',N,-ONE,
139 : $ NNZ,ROWPTR,COLIND,VAL,
140 : $ X,ONE,R)
141 : ENDIF
142 : *
143 : ITER = 0
144 :
145 : BNRM2 = DNRM2(N,B,1)
146 : IF(IMPR.GT.6) THEN
147 : WRITE(IOIMP,*) '||B||=',BNRM2
148 : ENDIF
149 :
150 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
151 : RESID = DNRM2(N,R,1) / BNRM2
152 : * RESID = DNRM2(N,R,1)
153 : IF(IMPR.GT.6) THEN
154 : WRITE(IOIMP,*) 'Résidu initial=',RESID
155 : ENDIF
156 :
157 : IF (RESID.LE.TOL) GOTO 30
158 : *
159 : 10 CONTINUE
160 : *
161 : * Construct the first column of V, and initialize S to the
162 : * elementary vector E1 scaled by RNORM.
163 : *
164 : C Preconditionner solve : Mv(1)=r
164 : CALL DCOPY(N,R,1,V(1,1),1)
165 : C s=||v(1)|| e(1)
166 : C v(1)=v(1) / ||v(1)||
167 : S(1) = DNRM2(N,V(1,1),1)
167 : CALL DSCAL(N,ONE/S(1),V(1,1),1)
168 : *
169 : DO 50 I = 1, RESTRT
170 : ITER = ITER + 1
171 : C av=A v(i)
171 : CALL DMOMV ('No transpose',N,ONE,
172 : $ NNZ,ROWPTR,COLIND,VAL,
173 : $ V(1,I),ZERO,AV)
174 : C Preconditionner solve : Mw=av
174 : CALL DCOPY(N,AV,1,W,1)
175 : *
176 : * Construct I-th column of H so that it is orthonormal to
177 : * the previous I-1 columns.
178 : *
178 : C CALL ORTHOH(I,N,RESTRT,WORK2(1,I+H-1),WORK(1,V),WORK(1,W))
179 : DO 60 K = 1, I
180 : H(K,I) = DDOT(N,W,1,V(1,K),1)
180 : CALL DAXPY(N,-H(K,I),V(1,K),1,W,1)
181 : 60 CONTINUE
182 : H(I+1,I)=DNRM2(N,W,1)
182 : CALL DCOPY(N,W,1,V(1,I+1),1)
182 : CALL DSCAL(N,ONE/H(I+1,I),V(1,I+1),1)
183 : *
184 : C IF (I.GT.1)
185 : *
186 : * Apply Givens rotations to the I-th column of H. This
187 : * effectively reduces the Hessenberg matrix to upper
188 : * triangular form during the RESTRT iterations.
189 : *
189 : C $ CALL APPGIV(I,RESTRT,WORK2(1,I+H-1),WORK2(1,GIV))
190 : DO 70 K = 1, I-1
191 : H1= (GCOS(K)*H(K,I)) + (GSIN(K)*H(K+1,I))
192 : H2=-(GSIN(K)*H(K,I)) + (GCOS(K)*H(K+1,I))
193 : H(K,I) =H1
194 : H(K+1,I)=H2
195 : 70 CONTINUE
196 : *
197 : * Approximate residual norm. Check tolerance. If okay, compute
198 : * final approximation vector X and quit.
199 : *
200 : RAY=SQRT((H(I,I)*H(I,I)) + (H(I+1,I)*H(I+1,I)))
201 : GCOS(I) = H(I,I) / RAY
202 : GSIN(I) = H(I+1,I) / RAY
203 : H(I,I)=RAY
204 : H(I+1,I)=ZERO
205 : S(I+1)=-GSIN(I)*S(I)
206 : S(I) = GCOS(I)*S(I)
207 : RESID=ABS(S(I+1)) / BNRM2
208 : LRES(ITER)=RESID
209 : IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN
210 : * Solve the system Hy=S
211 : * DO 501 J=I,1,-1
212 : * Y(J)=S(J) / H(J,J)
213 : * DO 502 K=J-1,1,-1
214 : * S(K)=S(K)-H(K,J)*Y(J)
215 : * 502 CONTINUE
216 : * 501 CONTINUE
216 : CALL DCOPY(I,S,1,Y,1)
216 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H,RESTRT+1,Y,1)
217 : * Updating solution
217 : C CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
218 : C $ WORK(1,S),WORK(1,V))
219 : DO 503 J=I,1,-1
219 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
220 : 503 CONTINUE
221 : IF (RESID.LE.TOL) THEN
222 : GOTO 30
223 : ENDIF
224 : IF (ITER.EQ.MAXIT) THEN
225 : GOTO 40
226 : ENDIF
227 : ENDIF
228 : 50 CONTINUE
229 : *
230 : * Compute current solution vector X.
231 : *
232 : * Solve the system Hy=S
233 : * DO 301 J=RESTRT,1,-1
234 : * Y(J)=S(J) / H(J,J)
235 : * DO 302 K=J-1,1,-1
236 : * S(K)=S(K)-H(K,J)*Y(J)
237 : * 302 CONTINUE
238 : * 301 CONTINUE
238 : CALL DCOPY(RESTRT,S,1,Y,1)
238 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',
239 : $ RESTRT,H,RESTRT+1,Y,1)
240 : * Updating solution
240 : C CALL UPDATE(RESTRT,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
241 : C $ WORK(1,S),WORK(1,V))
242 : DO 303 J=RESTRT,1,-1
242 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
243 : 303 CONTINUE
244 : *
245 : * Compute residual vector R, find norm,
246 : * then check for tolerance.
247 : *
247 : CALL DCOPY(N,B,1,R,1)
247 : CALL DMOMV ('No transpose',N,-ONE,
248 : $ NNZ,ROWPTR,COLIND,VAL,
249 : $ X,ONE,R)
250 : S(RESTRT+1) = DNRM2(N,R,1)
251 : RESID = S(RESTRT+1) / BNRM2
252 : LRES(ITER)=RESID
253 : IF (RESID.LE.TOL) GOTO 30
254 : *
255 : * Iteration fails.
256 : *
257 : IF (ITER.EQ.MAXIT) THEN
258 : GOTO 40
259 : ENDIF
260 : *
261 : * Do some more
262 : *
263 : GOTO 10
264 : *
265 : * Iteration successful
266 : *
267 : 30 CONTINUE
268 : IRET = 0
269 : RETURN
270 : *
271 : * Iteration fails
272 : *
273 : 40 CONTINUE
274 : IRET = 1
275 : RETURN
276 : *
277 : * End of GMR.
278 : *
279 : END
280 :
281 :
282 :
283 :
284 :
285 :
© Cast3M 2003 - All rights reserved.
Disclaimer