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