1 : C GMRDI SOURCE CHAT 05/01/13 00:18:41 5004
2 : SUBROUTINE GMRDI(N,NNZ,ROWPTR,COLIND,VAL,
3 : $ B,X,
4 : $ DINV,
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 : GMRDI
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 factorisation D-ILU de A.
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 PSDILU : Preconditionner solve
37 : C (Montée-Descente adaptée à D-ILU)
38 : C***********************************************************************
39 : C ENTREES : N, NNZ,
40 : C ROWPTR, COLIND, VAL,
41 : C B, DINV, RESTRT, IMPR
42 : C ENTREES/SORTIES : X,
43 : C R,W,AV,
44 : C Y,GCOS,GSIN,S,
45 : C V,H,
46 : C ITER, RESID
47 : C SORTIES : LRES, IRET
48 : C CODE RETOUR (IRET) : 0 si ok
49 : C <0 si problème (cf. BRTOL)
50 : C >0 si l'algorithme n'a pas convergé
51 : C N : nombre de degrés de liberté
52 : C NNZ : nombre de valeurs non nulles de la matrice Morse
53 : C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
54 : C et valeurs de la matrice Morse
55 : C B : vecteur second membre
56 : C DINV : vecteur contenant les inverses des pivots
57 : C de la factorisation D-ILU de A.
58 : C RESTRT : INTEGER
59 : C paramètre de redémarrage (restart) (i.e. m) pour
60 : C la méthode GMRES(m). (cf. prgmr.eso)
61 : C IMPR : niveau d'impression (0..4)
62 : C X : vecteur des inconnues.
63 : C En entrée, contient l'estimation x(0).
64 : C En sortie, la solution trouvée (si la méthode
65 : C a convergé).
66 : C R,W,AV : vecteurs de travail initialisés et détruits
67 : C Y,GCOS,GSIN,S dans prgmr.
68 : C V,H : tableaux de travail...
69 : C ITER : type INTEGER.
70 : C En entrée, contient le nombre maximum
71 : C d'itérations à effectuer.
72 : C En sortie, contient le nombre d'itérations
73 : C réellement effectuées.
74 : C RESID : type REAL*8.
75 : C En entrée, la valeur maximale du résidu à
76 : C convergence de l'algorithme mesurée comme suit :
77 : C norme[L2](b - A*x) / norme[L2]( b )
78 : C En sortie, la valeur effective de cette mesure.
79 : C LRES : pointeur sur segment MLREEL (include SMLREEL)
80 : C contient l'historique de la valeur de RESID en
81 : C fonction du nombre d'itérations effectuées.
82 : C
83 : C***********************************************************************
84 : C VERSION : v1, 16/06/98, version initiale
85 : C HISTORIQUE : v1, 16/06/98, création
86 : C HISTORIQUE :
87 : C HISTORIQUE :
88 : C***********************************************************************
89 : C Prière de PRENDRE LE TEMPS de compléter les commentaires
90 : C en cas de modification de ce sous-programme afin de faciliter
91 : C la maintenance !
92 : C***********************************************************************
93 : *
94 : * .. Includes .. a commenter si mise au point ok
95 : -INC CCOPTIO
96 : INTEGER IMPR
97 : * .. Scalar Arguments ..
98 : INTEGER N,NNZ
99 : INTEGER ITER,RESTRT,IRET
100 : REAL*8 RESID
101 : * ..
102 : * .. Array Arguments ..
103 : * .. Matrices stockées en Morse
104 : INTEGER ROWPTR(N+1)
105 : INTEGER COLIND(NNZ)
106 : REAL*8 VAL(NNZ)
107 : * .. Vecteurs
108 : REAL*8 B(N)
109 : REAL*8 X(N)
110 : REAL*8 DINV(N)
111 : REAL*8 LRES(ITER)
112 : * .. Tableaux de travail
113 : REAL*8 R(N),W(N),AV(N)
114 : REAL*8 Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT)
115 : REAL*8 S(RESTRT+1)
116 : REAL*8 V(N,RESTRT+1)
117 : REAL*8 H(RESTRT+1,RESTRT)
118 : *
119 : *
120 : * .. Variables locales
121 : * .. Parameters
122 : REAL*8 ZERO,ONE
123 : PARAMETER (ZERO=0.0D0, ONE=1.0D0)
124 : * ..
125 : INTEGER I,MAXIT,J,K
126 : REAL*8 BNRM2,TOL,DDOT,DNRM2
127 : REAL*8 H1,H2,RAY
128 : * ..
129 : * .. External subroutines and functions..
130 : EXTERNAL DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DTRSV
131 : EXTERNAL DMOMV,PSDILU
132 : * ..
133 : * .. Executable Statements ..
134 : *
135 : *
136 : IRET = 0
137 : MAXIT = ITER
138 : TOL = RESID
139 : IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmrdi.eso'
140 : *
141 : * Set initial residual.
142 : *
143 : C r=b
143 : CALL DCOPY(N,B,1,R,1)
144 : IF (DNRM2(N,X,1).NE.ZERO) THEN
145 : C r=b-Ax
145 : CALL DMOMV ('No transpose',N,-ONE,
146 : $ NNZ,ROWPTR,COLIND,VAL,
147 : $ X,ONE,R)
148 : ENDIF
149 : *
150 : ITER = 0
151 : *
152 : BNRM2 = DNRM2(N,B,1)
153 : IF(IMPR.GT.6) THEN
154 : WRITE(IOIMP,*) '||B||=',BNRM2
155 : ENDIF
156 :
157 : IF (BNRM2.EQ.ZERO) BNRM2 = ONE
158 : RESID = DNRM2(N,R,1) / BNRM2
159 : * RESID = DNRM2(N,R,1)
160 : IF(IMPR.GT.6) THEN
161 : WRITE(IOIMP,*) 'Résidu initial=',RESID
162 : ENDIF
163 :
164 : IF (RESID.LE.TOL) GOTO 30
165 : *
166 : 10 CONTINUE
167 : *
168 : * Construct the first column of V, and initialize S to the
169 : * elementary vector E1 scaled by RNORM.
170 : *
171 : C Preconditionner solve (D-ILU) : Mv(1)=r
171 : C CALL DCOPY(N,R,1,V(1,1),1)
171 : CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
172 : $ DINV,
173 : $ V(1,1),R)
174 :
175 : C s=||v(1)|| e(1)
176 : C v(1)=v(1) / ||v(1)||
177 : S(1) = DNRM2(N,V(1,1),1)
177 : CALL DSCAL(N,ONE/S(1),V(1,1),1)
178 : *
179 : DO 50 I = 1, RESTRT
180 : ITER = ITER + 1
181 : C av=A v(i)
181 : CALL DMOMV ('No transpose',N,ONE,
182 : $ NNZ,ROWPTR,COLIND,VAL,
183 : $ V(1,I),ZERO,AV)
184 : C Preconditionner solve (D-ILU) : Mw=av
184 : C CALL DCOPY(N,AV,1,W,1)
184 : CALL PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
185 : $ DINV,
186 : $ W,AV)
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 : *
210 : * Approximate residual norm. Check tolerance. If okay, compute
211 : * final approximation vector X and quit.
212 : *
213 : RAY=SQRT((H(I,I)*H(I,I)) + (H(I+1,I)*H(I+1,I)))
214 : GCOS(I) = H(I,I) / RAY
215 : GSIN(I) = H(I+1,I) / RAY
216 : H(I,I)=RAY
217 : H(I+1,I)=ZERO
218 : S(I+1)=-GSIN(I)*S(I)
219 : S(I) = GCOS(I)*S(I)
220 : RESID=ABS(S(I+1)) / BNRM2
221 : LRES(ITER)=RESID
222 : IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN
223 : * Solve the system Hy=S
224 : * DO 501 J=I,1,-1
225 : * Y(J)=S(J) / H(J,J)
226 : * DO 502 K=J-1,1,-1
227 : * S(K)=S(K)-H(K,J)*Y(J)
228 : * 502 CONTINUE
229 : * 501 CONTINUE
229 : CALL DCOPY(I,S,1,Y,1)
229 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H,RESTRT+1,Y,1)
230 : * Updating solution
230 : C CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
231 : C $ WORK(1,S),WORK(1,V))
232 : DO 503 J=I,1,-1
232 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
233 : 503 CONTINUE
234 : IF (RESID.LE.TOL) THEN
235 : GOTO 30
236 : ENDIF
237 : IF (ITER.EQ.MAXIT) THEN
238 : GOTO 40
239 : ENDIF
240 : ENDIF
241 : 50 CONTINUE
242 : *
243 : * Compute current solution vector X.
244 : *
245 : * Solve the system Hy=S
246 : * DO 301 J=RESTRT,1,-1
247 : * Y(J)=S(J) / H(J,J)
248 : * DO 302 K=J-1,1,-1
249 : * S(K)=S(K)-H(K,J)*Y(J)
250 : * 302 CONTINUE
251 : * 301 CONTINUE
251 : CALL DCOPY(RESTRT,S,1,Y,1)
251 : CALL DTRSV('UPPER','NOTRANS','NONUNIT',
252 : $ RESTRT,H,RESTRT+1,Y,1)
253 : * Updating solution
253 : C CALL UPDATE(RESTRT,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
254 : C $ WORK(1,S),WORK(1,V))
255 : DO 303 J=RESTRT,1,-1
255 : CALL DAXPY(N,Y(J),V(1,J),1,X,1)
256 : 303 CONTINUE
257 : *
258 : * Compute residual vector R, find norm,
259 : * then check for tolerance.
260 : *
260 : CALL DCOPY(N,B,1,R,1)
260 : CALL DMOMV ('No transpose',N,-ONE,
261 : $ NNZ,ROWPTR,COLIND,VAL,
262 : $ X,ONE,R)
263 : S(RESTRT+1) = DNRM2(N,R,1)
264 : RESID = S(RESTRT+1) / BNRM2
265 : LRES(ITER)=RESID
266 : IF (RESID.LE.TOL) GOTO 30
267 : *
268 : * Iteration fails.
269 : *
270 : IF (ITER.EQ.MAXIT) THEN
271 : GOTO 40
272 : ENDIF
273 : *
274 : * Do some more
275 : *
276 : GOTO 10
277 : *
278 : * Iteration successful
279 : *
280 : 30 CONTINUE
281 : IRET = 0
282 : RETURN
283 : *
284 : * Iteration fails
285 : *
286 : 40 CONTINUE
287 : IRET = 1
288 : RETURN
289 : *
290 : * End of GMRDI.
291 : *
292 : END
293 :
294 :
295 :
296 :
297 :
298 :
© Cast3M 2003 - All rights reserved.
Disclaimer