GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dspigm.f
Go to the documentation of this file.
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
subroutine datv(NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES, IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL, RPAR, IPAR)
Definition: datv.f:8
subroutine dhels(A, LDA, N, Q, B)
Definition: dhels.f:6
subroutine dheqr(A, LDA, N, Q, INFO, IJOB)
Definition: dheqr.f:6
subroutine dorth(VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
Definition: dorth.f:6
subroutine dspigm(NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL, MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V, HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS, RPAR, IPAR)
Definition: dspigm.f:9