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