|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DSPIGM (NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL, |
|
|
6 * MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V, |
|
|
7 * HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS, |
|
|
8 * RPAR, IPAR) |
|
|
9 C |
|
|
10 C***BEGIN PROLOGUE DSPIGM |
|
|
11 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
12 C***REVISION DATE 900926 (YYMMDD) |
|
|
13 C***REVISION DATE 940927 Removed MNEWT and added RHOK in call list. |
|
|
14 C |
|
|
15 C |
|
|
16 C----------------------------------------------------------------------- |
|
|
17 C***DESCRIPTION |
|
|
18 C |
|
|
19 C This routine solves the linear system A * Z = R using a scaled |
|
|
20 C preconditioned version of the generalized minimum residual method. |
|
|
21 C An initial guess of Z = 0 is assumed. |
|
|
22 C |
|
|
23 C On entry |
|
|
24 C |
|
|
25 C NEQ = Problem size, passed to PSOL. |
|
|
26 C |
|
|
27 C TN = Current Value of T. |
|
|
28 C |
|
|
29 C Y = Array Containing current dependent variable vector. |
|
|
30 C |
|
|
31 C YPRIME = Array Containing current first derivative of Y. |
|
|
32 C |
|
|
33 C SAVR = Array containing current value of G(T,Y,YPRIME). |
|
|
34 C |
|
|
35 C R = The right hand side of the system A*Z = R. |
|
|
36 C R is also used as work space when computing |
|
|
37 C the final approximation and will therefore be |
|
|
38 C destroyed. |
|
|
39 C (R is the same as V(*,MAXL+1) in the call to DSPIGM.) |
|
|
40 C |
|
|
41 C WGHT = The vector of length NEQ containing the nonzero |
|
|
42 C elements of the diagonal scaling matrix. |
|
|
43 C |
|
|
44 C MAXL = The maximum allowable order of the matrix H. |
|
|
45 C |
|
|
46 C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES. |
|
|
47 C |
|
|
48 C KMP = The number of previous vectors the new vector, VNEW, |
|
|
49 C must be made orthogonal to. (KMP .LE. MAXL.) |
|
|
50 C |
|
|
51 C EPLIN = Tolerance on residuals R-A*Z in weighted rms norm. |
|
|
52 C |
|
|
53 C CJ = Scalar proportional to current value of |
|
|
54 C 1/(step size H). |
|
|
55 C |
|
|
56 C WK = Real work array used by routine DATV and PSOL. |
|
|
57 C |
|
|
58 C DL = Real work array used for calculation of the residual |
|
|
59 C norm RHO when the method is incomplete (KMP.LT.MAXL) |
|
|
60 C and/or when using restarting. |
|
|
61 C |
|
|
62 C WP = Real work array used by preconditioner PSOL. |
|
|
63 C |
|
|
64 C IWP = Integer work array used by preconditioner PSOL. |
|
|
65 C |
|
|
66 C IRST = Method flag indicating if restarting is being |
|
|
67 C performed. IRST .GT. 0 means restarting is active, |
|
|
68 C while IRST = 0 means restarting is not being used. |
|
|
69 C |
|
|
70 C NRSTS = Counter for the number of restarts on the current |
|
|
71 C call to DSPIGM. If NRSTS .GT. 0, then the residual |
|
|
72 C R is already scaled, and so scaling of R is not |
|
|
73 C necessary. |
|
|
74 C |
|
|
75 C |
|
|
76 C On Return |
|
|
77 C |
|
|
78 C Z = The final computed approximation to the solution |
|
|
79 C of the system A*Z = R. |
|
|
80 C |
|
|
81 C LGMR = The number of iterations performed and |
|
|
82 C the current order of the upper Hessenberg |
|
|
83 C matrix HES. |
|
|
84 C |
|
|
85 C NRE = The number of calls to RES (i.e. DATV) |
|
|
86 C |
|
|
87 C NPSL = The number of calls to PSOL. |
|
|
88 C |
|
|
89 C V = The neq by (LGMR+1) array containing the LGMR |
|
|
90 C orthogonal vectors V(*,1) to V(*,LGMR). |
|
|
91 C |
|
|
92 C HES = The upper triangular factor of the QR decomposition |
|
|
93 C of the (LGMR+1) by LGMR upper Hessenberg matrix whose |
|
|
94 C entries are the scaled inner-products of A*V(*,I) |
|
|
95 C and V(*,K). |
|
|
96 C |
|
|
97 C Q = Real array of length 2*MAXL containing the components |
|
|
98 C of the givens rotations used in the QR decomposition |
|
|
99 C of HES. It is loaded in DHEQR and used in DHELS. |
|
|
100 C |
|
|
101 C IRES = Error flag from RES. |
|
|
102 C |
|
|
103 C DL = Scaled preconditioned residual, |
|
|
104 C (D-inverse)*(P-inverse)*(R-A*Z). Only loaded when |
|
|
105 C performing restarts of the Krylov iteration. |
|
|
106 C |
|
|
107 C RHOK = Weighted norm of final preconditioned residual. |
|
|
108 C |
|
|
109 C IFLAG = Integer error flag.. |
|
|
110 C 0 Means convergence in LGMR iterations, LGMR.LE.MAXL. |
|
|
111 C 1 Means the convergence test did not pass in MAXL |
|
|
112 C iterations, but the new residual norm (RHO) is |
|
|
113 C .LT. the old residual norm (RNRM), and so Z is |
|
|
114 C computed. |
|
|
115 C 2 Means the convergence test did not pass in MAXL |
|
|
116 C iterations, new residual norm (RHO) .GE. old residual |
|
|
117 C norm (RNRM), and the initial guess, Z = 0, is |
|
|
118 C returned. |
|
|
119 C 3 Means there was a recoverable error in PSOL |
|
|
120 C caused by the preconditioner being out of date. |
|
|
121 C -1 Means there was an unrecoverable error in PSOL. |
|
|
122 C |
|
|
123 C----------------------------------------------------------------------- |
|
|
124 C***ROUTINES CALLED |
|
|
125 C PSOL, DNRM2, DSCAL, DATV, DORTH, DHEQR, DCOPY, DHELS, DAXPY |
|
|
126 C |
|
|
127 C***END PROLOGUE DSPIGM |
|
|
128 C |
|
|
129 INTEGER NEQ,MAXL,MAXLP1,KMP,IRES,NRE,NPSL,LGMR,IWP, |
|
|
130 1 IFLAG,IRST,NRSTS,IPAR |
|
|
131 DOUBLE PRECISION TN,Y,YPRIME,SAVR,R,WGHT,EPLIN,CJ,Z,V,HES,Q,WP,WK, |
|
|
132 1 DL,RHOK,RPAR |
|
|
133 DIMENSION Y(*), YPRIME(*), SAVR(*), R(*), WGHT(*), Z(*), |
|
|
134 1 V(NEQ,*), HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*), |
|
|
135 2 RPAR(*), IPAR(*) |
|
|
136 INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1 |
|
|
137 DOUBLE PRECISION RNRM,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM |
|
|
138 EXTERNAL RES, PSOL |
|
|
139 C |
|
|
140 IER = 0 |
|
|
141 IFLAG = 0 |
|
|
142 LGMR = 0 |
|
|
143 NPSL = 0 |
|
|
144 NRE = 0 |
|
|
145 C----------------------------------------------------------------------- |
|
|
146 C The initial guess for Z is 0. The initial residual is therefore |
|
|
147 C the vector R. Initialize Z to 0. |
|
|
148 C----------------------------------------------------------------------- |
|
|
149 DO 10 I = 1,NEQ |
|
|
150 10 Z(I) = 0.0D0 |
|
|
151 C----------------------------------------------------------------------- |
|
|
152 C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0. |
|
|
153 C Form V(*,1), the scaled preconditioned right hand side. |
|
|
154 C----------------------------------------------------------------------- |
|
|
155 IF (NRSTS .EQ. 0) THEN |
|
|
156 CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, WK, CJ, WGHT, WP, IWP, |
|
|
157 1 R, EPLIN, IER, RPAR, IPAR) |
|
|
158 NPSL = 1 |
|
|
159 IF (IER .NE. 0) GO TO 300 |
|
|
160 DO 30 I = 1,NEQ |
|
|
161 30 V(I,1) = R(I)*WGHT(I) |
|
|
162 ELSE |
|
|
163 DO 35 I = 1,NEQ |
|
|
164 35 V(I,1) = R(I) |
|
|
165 ENDIF |
|
|
166 C----------------------------------------------------------------------- |
|
|
167 C Calculate norm of scaled vector V(*,1) and normalize it |
|
|
168 C If, however, the norm of V(*,1) (i.e. the norm of the preconditioned |
|
|
169 C residual) is .le. EPLIN, then return with Z=0. |
|
|
170 C----------------------------------------------------------------------- |
|
|
171 RNRM = DNRM2 (NEQ, V, 1) |
|
|
172 IF (RNRM .LE. EPLIN) THEN |
|
|
173 RHOK = RNRM |
|
|
174 RETURN |
|
|
175 ENDIF |
|
|
176 TEM = 1.0D0/RNRM |
|
|
177 CALL DSCAL (NEQ, TEM, V(1,1), 1) |
|
|
178 C----------------------------------------------------------------------- |
|
|
179 C Zero out the HES array. |
|
|
180 C----------------------------------------------------------------------- |
|
|
181 DO 65 J = 1,MAXL |
|
|
182 DO 60 I = 1,MAXLP1 |
|
|
183 60 HES(I,J) = 0.0D0 |
|
|
184 65 CONTINUE |
|
|
185 C----------------------------------------------------------------------- |
|
|
186 C Main loop to compute the vectors V(*,2) to V(*,MAXL). |
|
|
187 C The running product PROD is needed for the convergence test. |
|
|
188 C----------------------------------------------------------------------- |
|
|
189 PROD = 1.0D0 |
|
|
190 DO 90 LL = 1,MAXL |
|
|
191 LGMR = LL |
|
|
192 C----------------------------------------------------------------------- |
|
|
193 C Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is |
|
|
194 C the matrix A with scaling and inverse preconditioner factors applied. |
|
|
195 C Call routine DORTH to orthogonalize the new vector VNEW = V(*,LL+1). |
|
|
196 C call routine DHEQR to update the factors of HES. |
|
|
197 C----------------------------------------------------------------------- |
|
|
198 CALL DATV (NEQ, Y, TN, YPRIME, SAVR, V(1,LL), WGHT, Z, |
|
|
199 1 RES, IRES, PSOL, V(1,LL+1), WK, WP, IWP, CJ, EPLIN, |
|
|
200 1 IER, NRE, NPSL, RPAR, IPAR) |
|
|
201 IF (IRES .LT. 0) RETURN |
|
|
202 IF (IER .NE. 0) GO TO 300 |
|
|
203 CALL DORTH (V(1,LL+1), V, HES, NEQ, LL, MAXLP1, KMP, SNORMW) |
|
|
204 HES(LL+1,LL) = SNORMW |
|
|
205 CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL) |
|
|
206 IF (INFO .EQ. LL) GO TO 120 |
|
|
207 C----------------------------------------------------------------------- |
|
|
208 C Update RHO, the estimate of the norm of the residual R - A*ZL. |
|
|
209 C If KMP .LT. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not |
|
|
210 C necessarily orthogonal for LL .GT. KMP. The vector DL must then |
|
|
211 C be computed, and its norm used in the calculation of RHO. |
|
|
212 C----------------------------------------------------------------------- |
|
|
213 PROD = PROD*Q(2*LL) |
|
|
214 RHO = ABS(PROD*RNRM) |
|
|
215 IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN |
|
|
216 IF (LL .EQ. KMP+1) THEN |
|
|
217 CALL DCOPY (NEQ, V(1,1), 1, DL, 1) |
|
|
218 DO 75 I = 1,KMP |
|
|
219 IP1 = I + 1 |
|
|
220 I2 = I*2 |
|
|
221 S = Q(I2) |
|
|
222 C = Q(I2-1) |
|
|
223 DO 70 K = 1,NEQ |
|
|
224 70 DL(K) = S*DL(K) + C*V(K,IP1) |
|
|
225 75 CONTINUE |
|
|
226 ENDIF |
|
|
227 S = Q(2*LL) |
|
|
228 C = Q(2*LL-1)/SNORMW |
|
|
229 LLP1 = LL + 1 |
|
|
230 DO 80 K = 1,NEQ |
|
|
231 80 DL(K) = S*DL(K) + C*V(K,LLP1) |
|
|
232 DLNRM = DNRM2 (NEQ, DL, 1) |
|
|
233 RHO = RHO*DLNRM |
|
|
234 ENDIF |
|
|
235 C----------------------------------------------------------------------- |
|
|
236 C Test for convergence. If passed, compute approximation ZL. |
|
|
237 C If failed and LL .LT. MAXL, then continue iterating. |
|
|
238 C----------------------------------------------------------------------- |
|
|
239 IF (RHO .LE. EPLIN) GO TO 200 |
|
|
240 IF (LL .EQ. MAXL) GO TO 100 |
|
|
241 C----------------------------------------------------------------------- |
|
|
242 C Rescale so that the norm of V(1,LL+1) is one. |
|
|
243 C----------------------------------------------------------------------- |
|
|
244 TEM = 1.0D0/SNORMW |
|
|
245 CALL DSCAL (NEQ, TEM, V(1,LL+1), 1) |
|
|
246 90 CONTINUE |
|
|
247 100 CONTINUE |
|
|
248 IF (RHO .LT. RNRM) GO TO 150 |
|
|
249 120 CONTINUE |
|
|
250 IFLAG = 2 |
|
|
251 DO 130 I = 1,NEQ |
|
|
252 130 Z(I) = 0.D0 |
|
|
253 RETURN |
|
|
254 150 IFLAG = 1 |
|
|
255 C----------------------------------------------------------------------- |
|
|
256 C The tolerance was not met, but the residual norm was reduced. |
|
|
257 C If performing restarting (IRST .gt. 0) calculate the residual vector |
|
|
258 C RL and store it in the DL array. If the incomplete version is |
|
|
259 C being used (KMP .lt. MAXL) then DL has already been calculated. |
|
|
260 C----------------------------------------------------------------------- |
|
|
261 IF (IRST .GT. 0) THEN |
|
|
262 IF (KMP .EQ. MAXL) THEN |
|
|
263 C |
|
|
264 C Calculate DL from the V(I)'s. |
|
|
265 C |
|
|
266 CALL DCOPY (NEQ, V(1,1), 1, DL, 1) |
|
|
267 MAXLM1 = MAXL - 1 |
|
|
268 DO 175 I = 1,MAXLM1 |
|
|
269 IP1 = I + 1 |
|
|
270 I2 = I*2 |
|
|
271 S = Q(I2) |
|
|
272 C = Q(I2-1) |
|
|
273 DO 170 K = 1,NEQ |
|
|
274 170 DL(K) = S*DL(K) + C*V(K,IP1) |
|
|
275 175 CONTINUE |
|
|
276 S = Q(2*MAXL) |
|
|
277 C = Q(2*MAXL-1)/SNORMW |
|
|
278 DO 180 K = 1,NEQ |
|
|
279 180 DL(K) = S*DL(K) + C*V(K,MAXLP1) |
|
|
280 ENDIF |
|
|
281 C |
|
|
282 C Scale DL by RNRM*PROD to obtain the residual RL. |
|
|
283 C |
|
|
284 TEM = RNRM*PROD |
|
|
285 CALL DSCAL(NEQ, TEM, DL, 1) |
|
|
286 ENDIF |
|
|
287 C----------------------------------------------------------------------- |
|
|
288 C Compute the approximation ZL to the solution. |
|
|
289 C Since the vector Z was used as work space, and the initial guess |
|
|
290 C of the Newton correction is zero, Z must be reset to zero. |
|
|
291 C----------------------------------------------------------------------- |
|
|
292 200 CONTINUE |
|
|
293 LL = LGMR |
|
|
294 LLP1 = LL + 1 |
|
|
295 DO 210 K = 1,LLP1 |
|
|
296 210 R(K) = 0.0D0 |
|
|
297 R(1) = RNRM |
|
|
298 CALL DHELS (HES, MAXLP1, LL, Q, R) |
|
|
299 DO 220 K = 1,NEQ |
|
|
300 220 Z(K) = 0.0D0 |
|
|
301 DO 230 I = 1,LL |
|
|
302 CALL DAXPY (NEQ, R(I), V(1,I), 1, Z, 1) |
|
|
303 230 CONTINUE |
|
|
304 DO 240 I = 1,NEQ |
|
|
305 240 Z(I) = Z(I)/WGHT(I) |
|
|
306 C Load RHO into RHOK. |
|
|
307 RHOK = RHO |
|
|
308 RETURN |
|
|
309 C----------------------------------------------------------------------- |
|
|
310 C This block handles error returns forced by routine PSOL. |
|
|
311 C----------------------------------------------------------------------- |
|
|
312 300 CONTINUE |
|
|
313 IF (IER .LT. 0) IFLAG = -1 |
|
|
314 IF (IER .GT. 0) IFLAG = 3 |
|
|
315 C |
|
|
316 RETURN |
|
|
317 C |
|
|
318 C------END OF SUBROUTINE DSPIGM----------------------------------------- |
|
|
319 END |
