| blob | 11fc37a8f4aad49716f87f6e32ff9722420fe8f2 |
1 *DECK DLSODPK
8 C-----------------------------------------------------------------------
9 C This is the 18 November 2003 version of
10 C DLSODPK: Livermore Solver for Ordinary Differential equations,
11 C with Preconditioned Krylov iteration methods for the
12 C Newton correction linear systems.
13 C
14 C This version is in double precision.
15 C
16 C DLSODPK solves the initial value problem for stiff or nonstiff
17 C systems of first order ODEs,
18 C dy/dt = f(t,y) , or, in component form,
19 C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
20 C-----------------------------------------------------------------------
21 C Introduction.
22 C
23 C This is a modification of the DLSODE package which incorporates
24 C various preconditioned Krylov subspace iteration methods for the
25 C linear algebraic systems that arise in the case of stiff systems.
26 C
27 C The linear systems that must be solved have the form
28 C A * x = b , where A = identity - hl0 * (df/dy) .
29 C Here hl0 is a scalar, and df/dy is the Jacobian matrix of partial
30 C derivatives of f (NEQ by NEQ).
31 C
32 C The particular Krylov method is chosen by setting the second digit,
33 C MITER, in the method flag MF.
34 C Currently, the values of MITER have the following meanings:
35 C
36 C MITER = 1 means the preconditioned Scaled Incomplete
37 C Orthogonalization Method (SPIOM).
38 C
39 C 2 means an incomplete version of the Preconditioned Scaled
40 C Generalized Minimal Residual method (SPIGMR).
41 C This is the best choice in general.
42 C
43 C 3 means the Preconditioned Conjugate Gradient method (PCG).
44 C Recommended only when df/dy is symmetric or nearly so.
45 C
46 C 4 means the scaled Preconditioned Conjugate Gradient method
47 C (PCGS). Recommended only when D-inverse * df/dy * D is
48 C symmetric or nearly so, where D is the diagonal scaling
49 C matrix with elements 1/EWT(i) (see RTOL/ATOL description).
50 C
51 C 9 means that only a user-supplied matrix P (approximating A)
52 C will be used, with no Krylov iteration done. This option
53 C allows the user to provide the complete linear system
54 C solution algorithm, if desired.
55 C
56 C The user can apply preconditioning to the linear system A*x = b,
57 C by means of arbitrary matrices (the preconditioners).
58 C In the case of SPIOM and SPIGMR, one can apply left and right
59 C preconditioners P1 and P2, and the basic iterative method is then
60 C applied to the matrix (P1-inverse)*A*(P2-inverse) instead of to the
61 C matrix A. The product P1*P2 should be an approximation to matrix A
62 C such that linear systems with P1 or P2 are easier to solve than with
63 C A. Preconditioning from the left only or right only means using
64 C P2 = identity or P1 = identity, respectively.
65 C In the case of the PCG and PCGS methods, there is only one
66 C preconditioner matrix P (but it can be the product of more than one).
67 C It should approximate the matrix A but allow for relatively
68 C easy solution of linear systems with coefficient matrix P.
69 C For PCG, P should be positive definite symmetric, or nearly so,
70 C and for PCGS, the scaled preconditioner D-inverse * P * D
71 C should be symmetric or nearly so.
72 C If the Jacobian J = df/dy splits in a natural way into a sum
73 C J = J1 + J2, then one possible choice of preconditioners is
74 C P1 = identity - hl0 * J1 and P2 = identity - hl0 * J2
75 C provided each of these is easy to solve (or approximately solve).
76 C
77 C-----------------------------------------------------------------------
78 C References:
79 C 1. Peter N. Brown and Alan C. Hindmarsh, Reduced Storage Matrix
80 C Methods in Stiff ODE Systems, J. Appl. Math. & Comp., 31 (1989),
81 C pp. 40-91; also L.L.N.L. Report UCRL-95088, Rev. 1, June 1987.
82 C 2. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE
83 C Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.),
84 C North-Holland, Amsterdam, 1983, pp. 55-64.
85 C-----------------------------------------------------------------------
86 C Authors: Alan C. Hindmarsh and Peter N. Brown
87 C Center for Applied Scientific Computing, L-561
88 C Lawrence Livermore National Laboratory
89 C Livermore, CA 94551
90 C-----------------------------------------------------------------------
91 C Summary of Usage.
92 C
93 C Communication between the user and the DLSODPK package, for normal
94 C situations, is summarized here. This summary describes only a subset
95 C of the full set of options available. See the full description for
96 C details, including optional communication, nonstandard options,
97 C and instructions for special situations. See also the demonstration
98 C program distributed with this solver.
99 C
100 C A. First provide a subroutine of the form:
101 C SUBROUTINE F (NEQ, T, Y, YDOT)
102 C DOUBLE PRECISION T, Y(*), YDOT(*)
103 C which supplies the vector function f by loading YDOT(i) with f(i).
104 C
105 C B. Next determine (or guess) whether or not the problem is stiff.
106 C Stiffness occurs when the Jacobian matrix df/dy has an eigenvalue
107 C whose real part is negative and large in magnitude, compared to the
108 C reciprocal of the t span of interest. If the problem is nonstiff,
109 C use a method flag MF = 10. If it is stiff, MF should be between 21
110 C and 24, or possibly 29. MF = 22 is generally the best choice.
111 C Use 23 or 24 only if symmetry is present. Use MF = 29 if the
112 C complete linear system solution is to be provided by the user.
113 C The following four parameters must also be set.
114 C IWORK(1) = LWP = length of real array WP for preconditioning.
115 C IWORK(2) = LIWP = length of integer array IWP for preconditioning.
116 C IWORK(3) = JPRE = preconditioner type flag:
117 C = 0 for no preconditioning (P1 = P2 = P = identity)
118 C = 1 for left-only preconditioning (P2 = identity)
119 C = 2 for right-only preconditioning (P1 = identity)
120 C = 3 for two-sided preconditioning (and PCG or PCGS)
121 C IWORK(4) = JACFLG = flag for whether JAC is called.
122 C = 0 if JAC is not to be called,
123 C = 1 if JAC is to be called.
124 C Use JACFLG = 1 if JAC computes any nonconstant data for use in
125 C preconditioning, such as Jacobian elements.
126 C The arrays WP and IWP are work arrays under the user's control,
127 C for use in the routines that perform preconditioning operations.
128 C
129 C C. If the problem is stiff, you must supply two routines that deal
130 C with the preconditioning of the linear systems to be solved.
131 C These are as follows:
132 C
133 C SUBROUTINE JAC (F, NEQ, T, Y, YSV, REWT, FTY, V, HL0, WP,IWP, IER)
134 C DOUBLE PRECISION T, Y(*),YSV(*), REWT(*), FTY(*), V(*), HL0, WP(*)
135 C INTEGER IWP(*)
136 C This routine must evaluate and preprocess any parts of the
137 C Jacobian matrix df/dy involved in the preconditioners P1, P2, P.
138 C The Y and FTY arrays contain the current values of y and f(t,y),
139 C respectively, and YSV also contains the current value of y.
140 C The array V is work space of length NEQ.
141 C JAC must multiply all computed Jacobian elements by the scalar
142 C -HL0, add the identity matrix, and do any factorization
143 C operations called for, in preparation for solving linear systems
144 C with a coefficient matrix of P1, P2, or P. The matrix P1*P2 or P
145 C should be an approximation to identity - HL0 * (df/dy).
146 C JAC should return IER = 0 if successful, and IER .ne. 0 if not.
147 C (If IER .ne. 0, a smaller time step will be tried.)
148 C
149 C SUBROUTINE PSOL (NEQ, T, Y, FTY, WK, HL0, WP, IWP, B, LR, IER)
150 C DOUBLE PRECISION T, Y(*), FTY(*), WK(*), HL0, WP(*), B(*)
151 C INTEGER IWP(*)
152 C This routine must solve a linear system with B as right-hand
153 C side and one of the preconditioning matrices, P1, P2, or P, as
154 C coefficient matrix, and return the solution vector in B.
155 C LR is a flag concerning left vs right preconditioning, input
156 C to PSOL. PSOL is to use P1 if LR = 1 and P2 if LR = 2.
157 C In the case of the PCG or PCGS method, LR will be 3, and PSOL
158 C should solve the system P*x = B with the preconditioner matrix P.
159 C In the case MF = 29 (no Krylov iteration), LR will be 0,
160 C and PSOL is to return in B the desired approximate solution
161 C to A * x = B, where A = identity - HL0 * (df/dy).
162 C PSOL can use data generated in the JAC routine and stored in
163 C WP and IWP. WK is a work array of length NEQ.
164 C The argument HL0 is the current value of the scalar appearing
165 C in the linear system. If the old value, at the time of the last
166 C JAC call, is needed, it must have been saved by JAC in WP.
167 C On return, PSOL should set the error flag IER as follows:
168 C IER = 0 if PSOL was successful,
169 C IER .gt. 0 if a recoverable error occurred, meaning that the
170 C time step will be retried,
171 C IER .lt. 0 if an unrecoverable error occurred, meaning that the
172 C solver is to stop immediately.
173 C
174 C D. Write a main program which calls Subroutine DLSODPK once for
175 C each point at which answers are desired. This should also provide
176 C for possible use of logical unit 6 for output of error messages by
177 C DLSODPK. On the first call to DLSODPK, supply arguments as follows:
178 C F = name of subroutine for right-hand side vector f.
179 C This name must be declared External in calling program.
180 C NEQ = number of first order ODEs.
181 C Y = array of initial values, of length NEQ.
182 C T = the initial value of the independent variable.
183 C TOUT = first point where output is desired (.ne. T).
184 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
185 C RTOL = relative tolerance parameter (scalar).
186 C ATOL = absolute tolerance parameter (scalar or array).
187 C the estimated local error in y(i) will be controlled so as
188 C to be roughly less (in magnitude) than
189 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
190 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
191 C Thus the local error test passes if, in each component,
192 C either the absolute error is less than ATOL (or ATOL(i)),
193 C or the relative error is less than RTOL.
194 C Use RTOL = 0.0 for pure absolute error control, and
195 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
196 C control. Caution: Actual (global) errors may exceed these
197 C local tolerances, so choose them conservatively.
198 C ITASK = 1 for normal computation of output values of y at t = TOUT.
199 C ISTATE = integer flag (input and output). Set ISTATE = 1.
200 C IOPT = 0 to indicate no optional inputs used.
201 C RWORK = real work array of length at least:
202 C 20 + 16*NEQ for MF = 10,
203 C 45 + 17*NEQ + LWP for MF = 21,
204 C 61 + 17*NEQ + LWP for MF = 22,
205 C 20 + 15*NEQ + LWP for MF = 23 or 24,
206 C 20 + 12*NEQ + LWP for MF = 29.
207 C LRW = declared length of RWORK (in user's dimension).
208 C IWORK = integer work array of length at least:
209 C 30 for MF = 10,
210 C 35 + LIWP for MF = 21,
211 C 30 + LIWP for MF = 22, 23, 24, or 29.
212 C LIW = declared length of IWORK (in user's dimension).
213 C JAC,PSOL = names of subroutines for preconditioning.
214 C These names must be declared External in the calling program.
215 C MF = method flag. Standard values are:
216 C 10 for nonstiff (Adams) method.
217 C 21 for stiff (BDF) method, with preconditioned SIOM.
218 C 22 for stiff method, with preconditioned GMRES method.
219 C 23 for stiff method, with preconditioned CG method.
220 C 24 for stiff method, with scaled preconditioned CG method.
221 C 29 for stiff method, with user's PSOL routine only.
222 C Note that the main program must declare arrays Y, RWORK, IWORK,
223 C and possibly ATOL.
224 C
225 C E. The output from the first call (or any call) is:
226 C Y = array of computed values of y(t) vector.
227 C T = corresponding value of independent variable (normally TOUT).
228 C ISTATE = 2 if DLSODPK was successful, negative otherwise.
229 C -1 means excess work done on this call (perhaps wrong MF).
230 C -2 means excess accuracy requested (tolerances too small).
231 C -3 means illegal input detected (see printed message).
232 C -4 means repeated error test failures (check all inputs).
233 C -5 means repeated convergence failures (perhaps bad JAC
234 C or PSOL routine supplied or wrong choice of MF or
235 C tolerances, or this solver is inappropriate).
236 C -6 means error weight became zero during problem. (Solution
237 C component i vanished, and ATOL or ATOL(i) = 0.)
238 C -7 means an unrecoverable error occurred in PSOL.
239 C
240 C F. To continue the integration after a successful return, simply
241 C reset TOUT and call DLSODPK again. No other parameters need be reset.
242 C
243 C-----------------------------------------------------------------------
244 C-----------------------------------------------------------------------
245 C Full Description of User Interface to DLSODPK.
246 C
247 C The user interface to DLSODPK consists of the following parts.
248 C
249 C 1. The call sequence to Subroutine DLSODPK, which is a driver
250 C routine for the solver. This includes descriptions of both
251 C the call sequence arguments and of user-supplied routines.
252 C Following these descriptions is a description of
253 C optional inputs available through the call sequence, and then
254 C a description of optional outputs (in the work arrays).
255 C
256 C 2. Descriptions of other routines in the DLSODPK package that may be
257 C (optionally) called by the user. These provide the ability to
258 C alter error message handling, save and restore the internal
259 C Common, and obtain specified derivatives of the solution y(t).
260 C
261 C 3. Descriptions of Common blocks to be declared in overlay
262 C or similar environments, or to be saved when doing an interrupt
263 C of the problem and continued solution later.
264 C
265 C 4. Description of two routines in the DLSODPK package, either of
266 C which the user may replace with his/her own version, if desired.
267 C These relate to the measurement of errors.
268 C
269 C-----------------------------------------------------------------------
270 C Part 1. Call Sequence.
271 C
272 C The call sequence parameters used for input only are
273 C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, PSOL, MF,
274 C and those used for both input and output are
275 C Y, T, ISTATE.
276 C The work arrays RWORK and IWORK are also used for conditional and
277 C optional inputs and optional outputs. (The term output here refers
278 C to the return from Subroutine DLSODPK to the user's calling program.)
279 C
280 C The legality of input parameters will be thoroughly checked on the
281 C initial call for the problem, but not checked thereafter unless a
282 C change in input parameters is flagged by ISTATE = 3 on input.
283 C
284 C The descriptions of the call arguments are as follows.
285 C
286 C F = the name of the user-supplied subroutine defining the
287 C ODE system. The system must be put in the first-order
288 C form dy/dt = f(t,y), where f is a vector-valued function
289 C of the scalar t and the vector y. Subroutine F is to
290 C compute the function f. It is to have the form
291 C SUBROUTINE F (NEQ, T, Y, YDOT)
292 C DOUBLE PRECISION T, Y(*), YDOT(*)
293 C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
294 C is output. Y and YDOT are arrays of length NEQ.
295 C Subroutine F should not alter Y(1),...,Y(NEQ).
296 C F must be declared External in the calling program.
297 C
298 C Subroutine F may access user-defined quantities in
299 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
300 C (dimensioned in F) and/or Y has length exceeding NEQ(1).
301 C See the descriptions of NEQ and Y below.
302 C
303 C If quantities computed in the F routine are needed
304 C externally to DLSODPK, an extra call to F should be made
305 C for this purpose, for consistent and accurate results.
306 C If only the derivative dy/dt is needed, use DINTDY instead.
307 C
308 C NEQ = the size of the ODE system (number of first order
309 C ordinary differential equations). Used only for input.
310 C NEQ may be decreased, but not increased, during the problem.
311 C If NEQ is decreased (with ISTATE = 3 on input), the
312 C remaining components of Y should be left undisturbed, if
313 C these are to be accessed in the user-supplied subroutines.
314 C
315 C Normally, NEQ is a scalar, and it is generally referred to
316 C as a scalar in this user interface description. However,
317 C NEQ may be an array, with NEQ(1) set to the system size.
318 C (The DLSODPK package accesses only NEQ(1).) In either case,
319 C this parameter is passed as the NEQ argument in all calls
320 C to F, JAC, and PSOL. Hence, if it is an array, locations
321 C NEQ(2),... may be used to store other integer data and pass
322 C it to the user-supplied subroutines. Each such routine must
323 C include NEQ in a Dimension statement in that case.
324 C
325 C Y = a real array for the vector of dependent variables, of
326 C length NEQ or more. Used for both input and output on the
327 C first call (ISTATE = 1), and only for output on other calls.
328 C On the first call, Y must contain the vector of initial
329 C values. On output, Y contains the computed solution vector,
330 C evaluated at T. If desired, the Y array may be used
331 C for other purposes between calls to the solver.
332 C
333 C This array is passed as the Y argument in all calls to F,
334 C JAC, and PSOL. Hence its length may exceed NEQ, and locations
335 C Y(NEQ+1),... may be used to store other real data and
336 C pass it to the user-supplied subroutines. (The DLSODPK
337 C package accesses only Y(1),...,Y(NEQ).)
338 C
339 C T = the independent variable. On input, T is used only on the
340 C first call, as the initial point of the integration.
341 C On output, after each call, T is the value at which a
342 C computed solution y is evaluated (usually the same as TOUT).
343 C On an error return, T is the farthest point reached.
344 C
345 C TOUT = the next value of t at which a computed solution is desired.
346 C Used only for input.
347 C
348 C When starting the problem (ISTATE = 1), TOUT may be equal
349 C to T for one call, then should .ne. T for the next call.
350 C For the initial T, an input value of TOUT .ne. T is used
351 C in order to determine the direction of the integration
352 C (i.e. the algebraic sign of the step sizes) and the rough
353 C scale of the problem. Integration in either direction
354 C (forward or backward in t) is permitted.
355 C
356 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
357 C the first call (i.e. the first call with TOUT .ne. T).
358 C Otherwise, TOUT is required on every call.
359 C
360 C If ITASK = 1, 3, or 4, the values of TOUT need not be
361 C monotone, but a value of TOUT which backs up is limited
362 C to the current internal T interval, whose endpoints are
363 C TCUR - HU and TCUR (see optional outputs, below, for
364 C TCUR and HU).
365 C
366 C ITOL = an indicator for the type of error control. See
367 C description below under ATOL. Used only for input.
368 C
369 C RTOL = a relative error tolerance parameter, either a scalar or
370 C an array of length NEQ. See description below under ATOL.
371 C Input only.
372 C
373 C ATOL = an absolute error tolerance parameter, either a scalar or
374 C an array of length NEQ. Input only.
375 C
376 C The input parameters ITOL, RTOL, and ATOL determine
377 C the error control performed by the solver. The solver will
378 C control the vector E = (E(i)) of estimated local errors
379 C in y, according to an inequality of the form
380 C RMS-norm of ( E(i)/EWT(i) ) .le. 1,
381 C where EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
382 C and the RMS-norm (root-mean-square norm) here is
383 C RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
384 C is a vector of weights which must always be positive, and
385 C the values of RTOL and ATOL should all be non-negative.
386 C the following table gives the types (scalar/array) of
387 C RTOL and ATOL, and the corresponding form of EWT(i).
388 C
389 C ITOL RTOL ATOL EWT(i)
390 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
391 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
392 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
393 C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i)
394 C
395 C When either of these parameters is a scalar, it need not
396 C be dimensioned in the user's calling program.
397 C
398 C If none of the above choices (with ITOL, RTOL, and ATOL
399 C fixed throughout the problem) is suitable, more general
400 C error controls can be obtained by substituting
401 C user-supplied routines for the setting of EWT and/or for
402 C the norm calculation. See Part 4 below.
403 C
404 C If global errors are to be estimated by making a repeated
405 C run on the same problem with smaller tolerances, then all
406 C components of RTOL and ATOL (i.e. of EWT) should be scaled
407 C down uniformly.
408 C
409 C ITASK = an index specifying the task to be performed.
410 C Input only. ITASK has the following values and meanings.
411 C 1 means normal computation of output values of y(t) at
412 C t = TOUT (by overshooting and interpolating).
413 C 2 means take one step only and return.
414 C 3 means stop at the first internal mesh point at or
415 C beyond t = TOUT and return.
416 C 4 means normal computation of output values of y(t) at
417 C t = TOUT but without overshooting t = TCRIT.
418 C TCRIT must be input as RWORK(1). TCRIT may be equal to
419 C or beyond TOUT, but not behind it in the direction of
420 C integration. This option is useful if the problem
421 C has a singularity at or beyond t = TCRIT.
422 C 5 means take one step, without passing TCRIT, and return.
423 C TCRIT must be input as RWORK(1).
424 C
425 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
426 C (within roundoff), it will return T = TCRIT (exactly) to
427 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
428 C in which case answers at t = TOUT are returned first).
429 C
430 C ISTATE = an index used for input and output to specify the
431 C the state of the calculation.
432 C
433 C On input, the values of ISTATE are as follows.
434 C 1 means this is the first call for the problem
435 C (initializations will be done). See note below.
436 C 2 means this is not the first call, and the calculation
437 C is to continue normally, with no change in any input
438 C parameters except possibly TOUT and ITASK.
439 C (If ITOL, RTOL, and/or ATOL are changed between calls
440 C with ISTATE = 2, the new values will be used but not
441 C tested for legality.)
442 C 3 means this is not the first call, and the
443 C calculation is to continue normally, but with
444 C a change in input parameters other than
445 C TOUT and ITASK. Changes are allowed in
446 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF,
447 C and any of the optional inputs except H0.
448 C Note: A preliminary call with TOUT = T is not counted
449 C as a first call here, as no initialization or checking of
450 C input is done. (Such a call is sometimes useful for the
451 C purpose of outputting the initial conditions.)
452 C Thus the first call for which TOUT .ne. T requires
453 C ISTATE = 1 on input.
454 C
455 C On output, ISTATE has the following values and meanings.
456 C 1 means nothing was done; TOUT = T and ISTATE = 1 on input.
457 C 2 means the integration was performed successfully.
458 C -1 means an excessive amount of work (more than MXSTEP
459 C steps) was done on this call, before completing the
460 C requested task, but the integration was otherwise
461 C successful as far as T. (MXSTEP is an optional input
462 C and is normally 500.) To continue, the user may
463 C simply reset ISTATE to a value .gt. 1 and call again
464 C (the excess work step counter will be reset to 0).
465 C In addition, the user may increase MXSTEP to avoid
466 C this error return (see below on optional inputs).
467 C -2 means too much accuracy was requested for the precision
468 C of the machine being used. This was detected before
469 C completing the requested task, but the integration
470 C was successful as far as T. To continue, the tolerance
471 C parameters must be reset, and ISTATE must be set
472 C to 3. The optional output TOLSF may be used for this
473 C purpose. (Note: If this condition is detected before
474 C taking any steps, then an illegal input return
475 C (ISTATE = -3) occurs instead.)
476 C -3 means illegal input was detected, before taking any
477 C integration steps. See written message for details.
478 C Note: If the solver detects an infinite loop of calls
479 C to the solver with illegal input, it will cause
480 C the run to stop.
481 C -4 means there were repeated error test failures on
482 C one attempted step, before completing the requested
483 C task, but the integration was successful as far as T.
484 C The problem may have a singularity, or the input
485 C may be inappropriate.
486 C -5 means there were repeated convergence test failures on
487 C one attempted step, before completing the requested
488 C task, but the integration was successful as far as T.
489 C -6 means EWT(i) became zero for some i during the
490 C integration. Pure relative error control (ATOL(i)=0.0)
491 C was requested on a variable which has now vanished.
492 C The integration was successful as far as T.
493 C -7 means the PSOL routine returned an unrecoverable error
494 C flag (IER .lt. 0). The integration was successful as
495 C far as T.
496 C
497 C Note: since the normal output value of ISTATE is 2,
498 C it does not need to be reset for normal continuation.
499 C Also, since a negative input value of ISTATE will be
500 C regarded as illegal, a negative output value requires the
501 C user to change it, and possibly other inputs, before
502 C calling the solver again.
503 C
504 C IOPT = an integer flag to specify whether or not any optional
505 C inputs are being used on this call. Input only.
506 C The optional inputs are listed separately below.
507 C IOPT = 0 means no optional inputs are being used.
508 C Default values will be used in all cases.
509 C IOPT = 1 means one or more optional inputs are being used.
510 C
511 C RWORK = a real working array (double precision).
512 C The length of RWORK must be at least
513 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LENLS + LWP where
514 C NYH = the initial value of NEQ,
515 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
516 C smaller value is given as an optional input),
517 C LENLS = length of work space for linear system (Krylov)
518 C method, excluding preconditioning:
519 C LENLS = 0 if MITER = 0,
520 C LENLS = NEQ*(MAXL+3) + MAXL**2 if MITER = 1,
521 C LENLS = NEQ*(MAXL+3+MIN(1,MAXL-KMP))
522 C + (MAXL+3)*MAXL + 1 if MITER = 2,
523 C LENLS = 6*NEQ if MITER = 3 or 4,
524 C LENLS = 3*NEQ if MITER = 9.
525 C (See the MF description for METH and MITER, and the
526 C list of optional inputs for MAXL and KMP.)
527 C LWP = length of real user work space for preconditioning
528 C (see JAC/PSOL).
529 C Thus if default values are used and NEQ is constant,
530 C this length is:
531 C 20 + 16*NEQ for MF = 10,
532 C 45 + 24*NEQ + LWP FOR MF = 11,
533 C 61 + 24*NEQ + LWP FOR MF = 12,
534 C 20 + 22*NEQ + LWP FOR MF = 13 OR 14,
535 C 20 + 19*NEQ + LWP FOR MF = 19,
536 C 20 + 9*NEQ FOR MF = 20,
537 C 45 + 17*NEQ + LWP FOR MF = 21,
538 C 61 + 17*NEQ + LWP FOR MF = 22,
539 C 20 + 15*NEQ + LWP FOR MF = 23 OR 24,
540 C 20 + 12*NEQ + LWP for MF = 29.
541 C The first 20 words of RWORK are reserved for conditional
542 C and optional inputs and optional outputs.
543 C
544 C The following word in RWORK is a conditional input:
545 C RWORK(1) = TCRIT = critical value of t which the solver
546 C is not to overshoot. Required if ITASK is
547 C 4 or 5, and ignored otherwise. (See ITASK.)
548 C
549 C LRW = the length of the array RWORK, as declared by the user.
550 C (This will be checked by the solver.)
551 C
552 C IWORK = an integer work array. The length of IWORK must be at least
553 C 30 if MITER = 0 (MF = 10 or 20),
554 C 30 + MAXL + LIWP if MITER = 1 (MF = 11, 21),
555 C 30 + LIWP if MITER = 2, 3, 4, or 9.
556 C MAXL = 5 unless a different optional input value is given.
557 C LIWP = length of integer user work space for preconditioning
558 C (see conditional input list following).
559 C The first few words of IWORK are used for conditional and
560 C optional inputs and optional outputs.
561 C
562 C The following 4 words in IWORK are conditional inputs,
563 C required if MITER .ge. 1:
564 C IWORK(1) = LWP = length of real array WP for use in
565 C preconditioning (part of RWORK array).
566 C IWORK(2) = LIWP = length of integer array IWP for use in
567 C preconditioning (part of IWORK array).
568 C The arrays WP and IWP are work arrays under the
569 C user's control, for use in the routines that
570 C perform preconditioning operations (JAC and PSOL).
571 C IWORK(3) = JPRE = preconditioner type flag:
572 C = 0 for no preconditioning (P1 = P2 = P = identity)
573 C = 1 for left-only preconditioning (P2 = identity)
574 C = 2 for right-only preconditioning (P1 = identity)
575 C = 3 for two-sided preconditioning (and PCG or PCGS)
576 C IWORK(4) = JACFLG = flag for whether JAC is called.
577 C = 0 if JAC is not to be called,
578 C = 1 if JAC is to be called.
579 C Use JACFLG = 1 if JAC computes any nonconstant
580 C data needed in preconditioning operations,
581 C such as some of the Jacobian elements.
582 C
583 C LIW = the length of the array IWORK, as declared by the user.
584 C (This will be checked by the solver.)
585 C
586 C Note: The work arrays must not be altered between calls to DLSODPK
587 C for the same problem, except possibly for the conditional and
588 C optional inputs, and except for the last 3*NEQ words of RWORK.
589 C The latter space is used for internal scratch space, and so is
590 C available for use by the user outside DLSODPK between calls, if
591 C desired (but not for use by any of the user-supplied subroutines).
592 C
593 C JAC = the name of the user-supplied routine to compute any
594 C Jacobian elements (or approximations) involved in the
595 C matrix preconditioning operations (MITER .ge. 1).
596 C It is to have the form
597 C SUBROUTINE JAC (F, NEQ, T, Y, YSV, REWT, FTY, V,
598 C 1 HL0, WP, IWP, IER)
599 C DOUBLE PRECISION T, Y(*),YSV(*), REWT(*), FTY(*), V(*),
600 C 1 HL0, WP(*)
601 C INTEGER IWP(*)
602 C This routine must evaluate and preprocess any parts of the
603 C Jacobian matrix df/dy used in the preconditioners P1, P2, P.
604 C the Y and FTY arrays contain the current values of y and
605 C f(t,y), respectively, and YSV also contains the current
606 C value of y. The array V is work space of length
607 C NEQ for use by JAC. REWT is the array of reciprocal error
608 C weights (1/EWT). JAC must multiply all computed Jacobian
609 C elements by the scalar -HL0, add the identity matrix, and do
610 C any factorization operations called for, in preparation
611 C for solving linear systems with a coefficient matrix of
612 C P1, P2, or P. The matrix P1*P2 or P should be an
613 C approximation to identity - HL0 * (df/dy). JAC should
614 C return IER = 0 if successful, and IER .ne. 0 if not.
615 C (If IER .ne. 0, a smaller time step will be tried.)
616 C The arrays WP (of length LWP) and IWP (of length LIWP)
617 C are for use by JAC and PSOL for work space and for storage
618 C of data needed for the solution of the preconditioner
619 C linear systems. Their lengths and contents are under the
620 C user's control.
621 C The JAC routine may save relevant Jacobian elements (or
622 C approximations) used in the preconditioners, along with the
623 C value of HL0, and use these to reconstruct preconditioner
624 C matrices later without reevaluationg those elements.
625 C This may be cost-effective if JAC is called with HL0
626 C considerably different from its earlier value, indicating
627 C that a corrector convergence failure has occurred because
628 C of the change in HL0, not because of changes in the
629 C value of the Jacobian. In doing this, use the saved and
630 C current values of HL0 to decide whether to use saved
631 C or reevaluated elements.
632 C JAC may alter V, but may not alter Y, YSV, REWT, FTY, or HL0.
633 C JAC must be declared External in the calling program.
634 C Subroutine JAC may access user-defined quantities in
635 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
636 C (dimensioned in JAC) and/or Y has length exceeding NEQ(1).
637 C See the descriptions of NEQ and Y above.
638 C
639 C PSOL = the name of the user-supplied routine for the
640 C solution of preconditioner linear systems.
641 C It is to have the form
642 C SUBROUTINE PSOL (NEQ, T, Y, FTY, WK,HL0, WP,IWP, B, LR,IER)
643 C DOUBLE PRECISION T, Y(*), FTY(*), WK(*), HL0, WP(*), B(*)
644 C INTEGER IWP(*)
645 C This routine must solve a linear system with B as right-hand
646 C side and one of the preconditioning matrices, P1, P2, or P,
647 C as coefficient matrix, and return the solution vector in B.
648 C LR is a flag concerning left vs right preconditioning, input
649 C to PSOL. PSOL is to use P1 if LR = 1 and P2 if LR = 2.
650 C In the case of the PCG or PCGS method, LR will be 3, and PSOL
651 C should solve the system P*x = B with the preconditioner P.
652 C In the case MITER = 9 (no Krylov iteration), LR will be 0,
653 C and PSOL is to return in B the desired approximate solution
654 C to A * x = B, where A = identity - HL0 * (df/dy).
655 C PSOL can use data generated in the JAC routine and stored in
656 C WP and IWP.
657 C The Y and FTY arrays contain the current values of y and
658 C f(t,y), respectively. The array WK is work space of length
659 C NEQ for use by PSOL.
660 C The argument HL0 is the current value of the scalar appearing
661 C in the linear system. If the old value, as of the last
662 C JAC call, is needed, it must have been saved by JAC in WP.
663 C On return, PSOL should set the error flag IER as follows:
664 C IER = 0 if PSOL was successful,
665 C IER .gt. 0 on a recoverable error, meaning that the
666 C time step will be retried,
667 C IER .lt. 0 on an unrecoverable error, meaning that the
668 C solver is to stop immediately.
669 C PSOL may not alter Y, FTY, or HL0.
670 C PSOL must be declared External in the calling program.
671 C Subroutine PSOL may access user-defined quantities in
672 C NEQ(2),... and Y(NEQ(1)+1),... if NEQ is an array
673 C (dimensioned in PSOL) and/or Y has length exceeding NEQ(1).
674 C See the descriptions of NEQ and Y above.
675 C
676 C MF = the method flag. Used only for input. The legal values of
677 C MF are 10, 11, 12, 13, 14, 19, 20, 21, 22, 23, 24, and 29.
678 C MF has decimal digits METH and MITER: MF = 10*METH + MITER.
679 C METH indicates the basic linear multistep method:
680 C METH = 1 means the implicit Adams method.
681 C METH = 2 means the method based on Backward
682 C Differentiation Formulas (BDFs).
683 C MITER indicates the corrector iteration method:
684 C MITER = 0 means functional iteration (no linear system
685 C is involved).
686 C MITER = 1 means Newton iteration with Scaled Preconditioned
687 C Incomplete Orthogonalization Method (SPIOM)
688 C for the linear systems.
689 C MITER = 2 means Newton iteration with Scaled Preconditioned
690 C Generalized Minimal Residual method (SPIGMR)
691 C for the linear systems.
692 C MITER = 3 means Newton iteration with Preconditioned
693 C Conjugate Gradient method (PCG)
694 C for the linear systems.
695 C MITER = 4 means Newton iteration with scaled Preconditioned
696 C Conjugate Gradient method (PCGS)
697 C for the linear systems.
698 C MITER = 9 means Newton iteration with only the
699 C user-supplied PSOL routine called (no Krylov
700 C iteration) for the linear systems.
701 C JPRE is ignored, and PSOL is called with LR = 0.
702 C See comments in the introduction about the choice of MITER.
703 C If MITER .ge. 1, the user must supply routines JAC and PSOL
704 C (the names are arbitrary) as described above.
705 C For MITER = 0, dummy arguments can be used.
706 C-----------------------------------------------------------------------
707 C Optional Inputs.
708 C
709 C The following is a list of the optional inputs provided for in the
710 C call sequence. (See also Part 2.) For each such input variable,
711 C this table lists its name as used in this documentation, its
712 C location in the call sequence, its meaning, and the default value.
713 C The use of any of these inputs requires IOPT = 1, and in that
714 C case all of these inputs are examined. A value of zero for any
715 C of these optional inputs will cause the default value to be used.
716 C Thus to use a subset of the optional inputs, simply preload
717 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
718 C then set those of interest to nonzero values.
719 C
720 C Name Location Meaning and Default Value
721 C
722 C H0 RWORK(5) the step size to be attempted on the first step.
723 C The default value is determined by the solver.
724 C
725 C HMAX RWORK(6) the maximum absolute step size allowed.
726 C The default value is infinite.
727 C
728 C HMIN RWORK(7) the minimum absolute step size allowed.
729 C The default value is 0. (This lower bound is not
730 C enforced on the final step before reaching TCRIT
731 C when ITASK = 4 or 5.)
732 C
733 C DELT RWORK(8) convergence test constant in Krylov iteration
734 C algorithm. The default is .05.
735 C
736 C MAXORD IWORK(5) the maximum order to be allowed. The default
737 C value is 12 if METH = 1, and 5 if METH = 2.
738 C If MAXORD exceeds the default value, it will
739 C be reduced to the default value.
740 C If MAXORD is changed during the problem, it may
741 C cause the current order to be reduced.
742 C
743 C MXSTEP IWORK(6) maximum number of (internally defined) steps
744 C allowed during one call to the solver.
745 C The default value is 500.
746 C
747 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
748 C warning that T + H = T on a step (H = step size).
749 C This must be positive to result in a non-default
750 C value. The default value is 10.
751 C
752 C MAXL IWORK(8) maximum number of iterations in the SPIOM, SPIGMR,
753 C PCG, or PCGS algorithm (.le. NEQ).
754 C The default is MAXL = MIN(5,NEQ).
755 C
756 C KMP IWORK(9) number of vectors on which orthogonalization
757 C is done in SPIOM or SPIGMR algorithm (.le. MAXL).
758 C The default is KMP = MAXL.
759 C Note: When KMP .lt. MAXL and MF = 22, the length
760 C of RWORK must be defined accordingly. See
761 C the definition of RWORK above.
762 C-----------------------------------------------------------------------
763 C Optional Outputs.
764 C
765 C As optional additional output from DLSODPK, the variables listed
766 C below are quantities related to the performance of DLSODPK
767 C which are available to the user. These are communicated by way of
768 C the work arrays, but also have internal mnemonic names as shown.
769 C Except where stated otherwise, all of these outputs are defined
770 C on any successful return from DLSODPK, and on any return with
771 C ISTATE = -1, -2, -4, -5, -6, or -7. On an illegal input return
772 C (ISTATE = -3), they will be unchanged from their existing values
773 C (if any), except possibly for TOLSF, LENRW, and LENIW.
774 C On any error return, outputs relevant to the error will be defined,
775 C as noted below.
776 C
777 C Name Location Meaning
778 C
779 C HU RWORK(11) the step size in t last used (successfully).
780 C
781 C HCUR RWORK(12) the step size to be attempted on the next step.
782 C
783 C TCUR RWORK(13) the current value of the independent variable
784 C which the solver has actually reached, i.e. the
785 C current internal mesh point in t. On output, TCUR
786 C will always be at least as far as the argument
787 C T, but may be farther (if interpolation was done).
788 C
789 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
790 C computed when a request for too much accuracy was
791 C detected (ISTATE = -3 if detected at the start of
792 C the problem, ISTATE = -2 otherwise). If ITOL is
793 C left unaltered but RTOL and ATOL are uniformly
794 C scaled up by a factor of TOLSF for the next call,
795 C then the solver is deemed likely to succeed.
796 C (The user may also ignore TOLSF and alter the
797 C tolerance parameters in any other way appropriate.)
798 C
799 C NST IWORK(11) the number of steps taken for the problem so far.
800 C
801 C NFE IWORK(12) the number of f evaluations for the problem so far.
802 C
803 C NPE IWORK(13) the number of calls to JAC so far (for Jacobian
804 C evaluation associated with preconditioning).
805 C
806 C NQU IWORK(14) the method order last used (successfully).
807 C
808 C NQCUR IWORK(15) the order to be attempted on the next step.
809 C
810 C IMXER IWORK(16) the index of the component of largest magnitude in
811 C the weighted local error vector ( E(i)/EWT(i) ),
812 C on an error return with ISTATE = -4 or -5.
813 C
814 C LENRW IWORK(17) the length of RWORK actually required.
815 C This is defined on normal returns and on an illegal
816 C input return for insufficient storage.
817 C
818 C LENIW IWORK(18) the length of IWORK actually required.
819 C This is defined on normal returns and on an illegal
820 C input return for insufficient storage.
821 C
822 C NNI IWORK(19) number of nonlinear iterations so far (each of
823 C which calls an iterative linear solver).
824 C
825 C NLI IWORK(20) number of linear iterations so far.
826 C Note: A measure of the success of algorithm is
827 C the average number of linear iterations per
828 C nonlinear iteration, given by NLI/NNI.
829 C If this is close to MAXL, MAXL may be too small.
830 C
831 C NPS IWORK(21) number of preconditioning solve operations
832 C (PSOL calls) so far.
833 C
834 C NCFN IWORK(22) number of convergence failures of the nonlinear
835 C (Newton) iteration so far.
836 C Note: A measure of success is the overall
837 C rate of nonlinear convergence failures, NCFN/NST.
838 C
839 C NCFL IWORK(23) number of convergence failures of the linear
840 C iteration so far.
841 C Note: A measure of success is the overall
842 C rate of linear convergence failures, NCFL/NNI.
843 C
844 C The following two arrays are segments of the RWORK array which
845 C may also be of interest to the user as optional outputs.
846 C For each array, the table below gives its internal name,
847 C its base address in RWORK, and its description.
848 C
849 C Name Base Address Description
850 C
851 C YH 21 the Nordsieck history array, of size NYH by
852 C (NQCUR + 1), where NYH is the initial value
853 C of NEQ. For j = 0,1,...,NQCUR, column j+1
854 C of YH contains HCUR**j/factorial(j) times
855 C the j-th derivative of the interpolating
856 C polynomial currently representing the solution,
857 C evaluated at t = TCUR.
858 C
859 C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
860 C corrections on each step, scaled on output
861 C to represent the estimated local error in y
862 C on the last step. This is the vector E in
863 C the description of the error control. It is
864 C defined only on a successful return from
865 C DLSODPK.
866 C
867 C-----------------------------------------------------------------------
868 C Part 2. Other Routines Callable.
869 C
870 C The following are optional calls which the user may make to
871 C gain additional capabilities in conjunction with DLSODPK.
872 C (The routines XSETUN and XSETF are designed to conform to the
873 C SLATEC error handling package.)
874 C
875 C Form of Call Function
876 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
877 C output of messages from DLSODPK, if
878 C the default is not desired.
879 C The default value of lun is 6.
880 C
881 C CALL XSETF(MFLAG) Set a flag to control the printing of
882 C messages by DLSODPK.
883 C MFLAG = 0 means do not print. (Danger:
884 C This risks losing valuable information.)
885 C MFLAG = 1 means print (the default).
886 C
887 C Either of the above calls may be made at
888 C any time and will take effect immediately.
889 C
890 C CALL DSRCPK(RSAV,ISAV,JOB) saves and restores the contents of
891 C the internal Common blocks used by
892 C DLSODPK (see Part 3 below).
893 C RSAV must be a real array of length 222
894 C or more, and ISAV must be an integer
895 C array of length 50 or more.
896 C JOB=1 means save Common into RSAV/ISAV.
897 C JOB=2 means restore Common from RSAV/ISAV.
898 C DSRCPK is useful if one is
899 C interrupting a run and restarting
900 C later, or alternating between two or
901 C more problems solved with DLSODPK.
902 C
903 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
904 C (See below) orders, at a specified point t, if
905 C desired. It may be called only after
906 C a successful return from DLSODPK.
907 C
908 C The detailed instructions for using DINTDY are as follows.
909 C The form of the call is:
910 C
911 C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
912 C
913 C The input parameters are:
914 C
915 C T = value of independent variable where answers are desired
916 C (normally the same as the T last returned by DLSODPK).
917 C for valid results, T must lie between TCUR - HU and TCUR.
918 C (See optional outputs for TCUR and HU.)
919 C K = integer order of the derivative desired. K must satisfy
920 C 0 .le. K .le. NQCUR, where NQCUR is the current order
921 C (see optional outputs). The capability corresponding
922 C to K = 0, i.e. computing y(T), is already provided
923 C by DLSODPK directly. Since NQCUR .ge. 1, the first
924 C derivative dy/dt is always available with DINTDY.
925 C RWORK(21) = the base address of the history array YH.
926 C NYH = column length of YH, equal to the initial value of NEQ.
927 C
928 C The output parameters are:
929 C
930 C DKY = a real array of length NEQ containing the computed value
931 C of the K-th derivative of y(t).
932 C IFLAG = integer flag, returned as 0 if K and T were legal,
933 C -1 if K was illegal, and -2 if T was illegal.
934 C On an error return, a message is also written.
935 C-----------------------------------------------------------------------
936 C Part 3. Common Blocks.
937 C
938 C If DLSODPK is to be used in an overlay situation, the user
939 C must declare, in the primary overlay, the variables in:
940 C (1) the call sequence to DLSODPK, and
941 C (2) the two internal Common blocks
942 C /DLS001/ of length 255 (218 double precision words
943 C followed by 37 integer words),
944 C /DLPK01/ of length 17 (4 double precision words
945 C followed by 13 integer words).
946 C
947 C If DLSODPK is used on a system in which the contents of internal
948 C Common blocks are not preserved between calls, the user should
949 C declare the above Common blocks in the calling program to insure
950 C that their contents are preserved.
951 C
952 C If the solution of a given problem by DLSODPK is to be interrupted
953 C and then later continued, such as when restarting an interrupted run
954 C or alternating between two or more problems, the user should save,
955 C following the return from the last DLSODPK call prior to the
956 C interruption, the contents of the call sequence variables and the
957 C internal Common blocks, and later restore these values before the
958 C next DLSODPK call for that problem. To save and restore the Common
959 C blocks, use Subroutine DSRCPK (see Part 2 above).
960 C
961 C-----------------------------------------------------------------------
962 C Part 4. Optionally Replaceable Solver Routines.
963 C
964 C below are descriptions of two routines in the DLSODPK package which
965 C relate to the measurement of errors. Either routine can be
966 C replaced by a user-supplied version, if desired. However, since such
967 C a replacement may have a major impact on performance, it should be
968 C done only when absolutely necessary, and only with great caution.
969 C (Note: The means by which the package version of a routine is
970 C superseded by the user's version may be system-dependent.)
971 C
972 C (a) DEWSET.
973 C The following subroutine is called just before each internal
974 C integration step, and sets the array of error weights, EWT, as
975 C described under ITOL/RTOL/ATOL above:
976 C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
977 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODPK call sequence,
978 C YCUR contains the current dependent variable vector, and
979 C EWT is the array of weights set by DEWSET.
980 C
981 C If the user supplies this subroutine, it must return in EWT(i)
982 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
983 C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
984 C routine (see below), and also used by DLSODPK in the computation
985 C of the optional output IMXER, the diagonal Jacobian approximation,
986 C and the increments for difference quotient Jacobians.
987 C
988 C In the user-supplied version of DEWSET, it may be desirable to use
989 C the current values of derivatives of y. Derivatives up to order NQ
990 C are available from the history array YH, described above under
991 C optional outputs. In DEWSET, YH is identical to the YCUR array,
992 C extended to NQ + 1 columns with a column length of NYH and scale
993 C factors of H**j/factorial(j). On the first call for the problem,
994 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
995 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
996 C can be obtained by including in DEWSET the statements:
997 C DOUBLE PRECISION RLS
998 C COMMON /DLS001/ RLS(218),ILS(37)
999 C NQ = ILS(33)
1000 C NST = ILS(34)
1001 C H = RLS(212)
1002 C Thus, for example, the current value of dy/dt can be obtained as
1003 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
1004 C unnecessary when NST = 0).
1005 C
1006 C (b) DVNORM.
1007 C The following is a real function routine which computes the weighted
1008 C root-mean-square norm of a vector v:
1009 C D = DVNORM (N, V, W)
1010 C where:
1011 C N = the length of the vector,
1012 C V = real array of length N containing the vector,
1013 C W = real array of length N containing weights,
1014 C D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).
1015 C DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
1016 C EWT is as set by Subroutine DEWSET.
1017 C
1018 C If the user supplies this function, it should return a non-negative
1019 C value of DVNORM suitable for use in the error control in DLSODPK.
1020 C None of the arguments should be altered by DVNORM.
1021 C For example, a user-supplied DVNORM routine might:
1022 C -substitute a max-norm of (V(i)*W(i)) for the RMS-norm, or
1023 C -ignore some components of V in the norm, with the effect of
1024 C suppressing the error control on those components of y.
1025 C-----------------------------------------------------------------------
1026 C
1027 C***REVISION HISTORY (YYYYMMDD)
1028 C 19860901 DATE WRITTEN
1029 C 19861010 Numerous minor revisions to SPIOM and SPGMR routines;
1030 C minor corrections to prologues and comments.
1031 C 19870114 Changed name SPGMR to SPIGMR; revised residual norm
1032 C calculation in SPIGMR (for incomplete case);
1033 C revised error return logic in SPIGMR;
1034 C 19870330 Major update: corrected comments throughout;
1035 C removed TRET from Common; rewrote EWSET with 4 loops;
1036 C fixed t test in INTDY; added Cray directives in STODPK;
1037 C in STODPK, fixed DELP init. and logic around PJAC call;
1038 C combined routines to save/restore Common;
1039 C passed LEVEL = 0 in error message calls (except run abort).
1040 C 19871130 Added option MITER = 9; shortened WM array by 2;
1041 C revised early return from SPIOM and SPIGMR;
1042 C replaced copy loops with SCOPY/DCOPY calls;
1043 C minor corrections/revisions to SOLPK, SPIGMR, ATV, ATP;
1044 C corrections to main prologue and internal comments.
1045 C 19880304 Corrections to type declarations in SOLPK, SPIOM, USOL.
1046 C 19891025 Added ISTATE = -7 return; minor revisions to USOL;
1047 C added initialization of JACFLG in main driver;
1048 C removed YH and NYH from PKSET call list;
1049 C minor revisions to SPIOM and SPIGMR;
1050 C corrections to main prologue and internal comments.
1051 C 19900803 Added YSV to JAC call list; minor comment corrections.
1052 C 20010425 Major update: convert source lines to upper case;
1053 C added *DECK lines; changed from 1 to * in dummy dimensions;
1054 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
1055 C renamed routines for uniqueness across single/double prec.;
1056 C converted intrinsic names to generic form;
1057 C removed ILLIN and NTREP (data loaded) from Common;
1058 C removed all 'own' variables from Common;
1059 C changed error messages to quoted strings;
1060 C replaced XERRWV/XERRWD with 1993 revised version;
1061 C converted prologues, comments, error messages to mixed case;
1062 C numerous corrections to prologues and internal comments.
1063 C 20010507 Converted single precision source to double precision.
1064 C 20020502 Corrected declarations in descriptions of user routines.
1065 C 20030603 Corrected duplicate type declaration for DUMACH.
1066 C 20031105 Restored 'own' variables to Common blocks, to enable
1067 C interrupt/restart feature.
1068 C 20031112 Added SAVE statements for data-loaded constants.
1069 C 20031117 Changed internal name NPE to NJE.
1070 C
1071 C-----------------------------------------------------------------------
1072 C Other routines in the DLSODPK package.
1073 C
1074 C In addition to Subroutine DLSODPK, the DLSODPK package includes the
1075 C following subroutines and function routines:
1076 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1077 C DEWSET sets the error weight vector EWT before each step.
1078 C DVNORM computes the weighted RMS-norm of a vector.
1079 C DSTODPK is the core integrator, which does one step of the
1080 C integration and the associated error control.
1081 C DCFODE sets all method coefficients and test constants.
1082 C DPKSET interfaces between DSTODPK and the JAC routine.
1083 C DSOLPK manages solution of linear system in Newton iteration.
1084 C DSPIOM performs the SPIOM algorithm.
1085 C DATV computes a scaled, preconditioned product (I-hl0*J)*v.
1086 C DORTHOG orthogonalizes a vector against previous basis vectors.
1087 C DHEFA generates an LU factorization of a Hessenberg matrix.
1088 C DHESL solves a Hessenberg square linear system.
1089 C DSPIGMR performs the SPIGMR algorithm.
1090 C DHEQR generates a QR factorization of a Hessenberg matrix.
1091 C DHELS finds the least squares solution of a Hessenberg system.
1092 C DPCG performs Preconditioned Conjugate Gradient algorithm (PCG).
1093 C DPCGS performs the PCGS algorithm.
1094 C DATP computes the product A*p, where A = I - hl0*df/dy.
1095 C DUSOL interfaces to the user's PSOL routine (MITER = 9).
1096 C DSRCPK is a user-callable routine to save and restore
1097 C the contents of the internal Common blocks.
1098 C DAXPY, DCOPY, DDOT, DNRM2, and DSCAL are basic linear
1099 C algebra modules (from the BLAS collection).
1100 C DUMACH computes the unit roundoff in a machine-independent manner.
1101 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
1102 C error messages and warnings. XERRWD is machine-dependent.
1103 C Note: DVNORM, DDOT, DNRM2, DUMACH, IXSAV, and IUMACH are function
1104 C routines. All the others are subroutines.
1105 C
1106 C-----------------------------------------------------------------------
1125 CHARACTER*60 MSG
1127 C-----------------------------------------------------------------------
1128 C The following two internal Common blocks contain
1129 C (a) variables which are local to any subroutine but whose values must
1130 C be preserved between calls to the routine ("own" variables), and
1131 C (b) variables which are communicated between subroutines.
1132 C The block DLS001 is declared in subroutines DLSODPK, DINTDY, DSTODPK,
1133 C DSOLPK, and DATV.
1134 C The block DLPK01 is declared in subroutines DLSODPK, DSTODPK, DPKSET,
1135 C and DSOLPK.
1136 C Groups of variables are replaced by dummy arrays in the Common
1137 C declarations in routines where those variables are not used.
1138 C-----------------------------------------------------------------------
1145 C
1149 C
1151 C-----------------------------------------------------------------------
1152 C Block A.
1153 C This code block is executed on every call.
1154 C It tests ISTATE and ITASK for legality and branches appropriately.
1155 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1156 C not yet been done, an error return occurs.
1157 C If ISTATE = 1 and TOUT = T, return immediately.
1158 C-----------------------------------------------------------------------
1167 C-----------------------------------------------------------------------
1168 C Block B.
1169 C The next code block is executed for the initial call (ISTATE = 1),
1170 C or for a continuation call with parameter changes (ISTATE = 3).
1171 C It contains checking of all inputs and various initializations.
1172 C
1173 C First check legality of the non-optional inputs NEQ, ITOL, IOPT, MF.
1174 C-----------------------------------------------------------------------
1189 C Next process and check the optional inputs. --------------------------
1192 MXSTEP = MXSTP0
1193 MXHNIL = MXHNL0
1198 KMP = MAXL
1227 C-----------------------------------------------------------------------
1228 C Set work array pointers and check lengths LRW and LIW.
1229 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1230 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1231 C RWORK segments (in order) are denoted YH, WM, EWT, SAVF, SAVX, ACOR.
1232 C-----------------------------------------------------------------------
1263 C Check RTOL and ATOL for legality. ------------------------------------
1272 C Load SQRT(N) and its reciprocal in Common. ---------------------------
1276 C If ISTATE = 3, set flag to signal parameter changes to DSTODPK. ------
1279 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1284 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1291 C-----------------------------------------------------------------------
1292 C Block C.
1293 C The next block is for the initial call only (ISTATE = 1).
1294 C It contains all remaining initializations, the initial call to F,
1295 C and the calculation of the initial step size.
1296 C The error weights in EWT are inverted after being loaded.
1297 C-----------------------------------------------------------------------
1299 TN = T
1326 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1330 C Load the initial value vector in YH. ---------------------------------
1333 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1340 C-----------------------------------------------------------------------
1341 C The coding below computes the step size, H0, to be attempted on the
1342 C first step, unless the user has supplied a value for this.
1343 C First check that TOUT - T differs significantly from zero.
1344 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1345 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1346 C so as to be between 100*UROUND and 1.0E-3.
1347 C Then the computed value H0 is given by..
1348 C NEQ
1349 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * Sum ( f(i)/ywt(i) )**2 )
1350 C 1
1351 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1352 C f(i) = i-th component of initial value of f,
1353 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1354 C The sign of H0 is inferred from the initial values of TOUT and T.
1355 C-----------------------------------------------------------------------
1378 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1381 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1382 H = H0
1386 C-----------------------------------------------------------------------
1387 C Block D.
1388 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1389 C and is to check stop conditions before taking a step.
1390 C-----------------------------------------------------------------------
1392 NLI0 = NLI
1393 NNI0 = NNI
1394 NCFN0 = NCFN
1395 NCFL0 = NCFL
1401 T = TOUT
1413 T = TOUT
1424 C-----------------------------------------------------------------------
1425 C Block E.
1426 C The next block is normally executed for all calls and contains
1427 C the call to the one-step core integrator DSTODPK.
1428 C
1429 C This is a looping point for the integration steps.
1430 C
1431 C First check for too many steps being taken,
1432 C Check for poor Newton/Krylov method performance, update EWT (if not
1433 C at start of problem), check for too much accuracy being requested,
1434 C and check for H below the roundoff level in T.
1435 C-----------------------------------------------------------------------
1454 ENDIF
1458 ENDIF
1462 ENDIF
1466 ENDIF
1470 ENDIF
1474 ENDIF
1500 C-----------------------------------------------------------------------
1501 C CALL DSTODPK(NEQ,Y,YH,NYH,YH,EWT,SAVF,SAVX,ACOR,WM,IWM,F,JAC,PSOL)
1502 C-----------------------------------------------------------------------
1508 C-----------------------------------------------------------------------
1509 C Block F.
1510 C The following block handles the case of a successful return from the
1511 C core integrator (KFLAG = 0). Test for stop conditions.
1512 C-----------------------------------------------------------------------
1515 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1518 T = TOUT
1520 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1523 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1526 T = TOUT
1536 C ITASK = 5. see if TCRIT was reached and jump to exit. ---------------
1539 C-----------------------------------------------------------------------
1540 C Block G.
1541 C The following block handles all successful returns from DLSODPK.
1542 C If ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1543 C ISTATE is set to 2, and the optional outputs are loaded into the
1544 C work arrays before returning.
1545 C-----------------------------------------------------------------------
1548 T = TN
1565 RETURN
1566 C-----------------------------------------------------------------------
1567 C Block H.
1568 C The following block handles all unsuccessful returns other than
1569 C those for illegal input. First the error message routine is called.
1570 C If there was an error test or convergence test failure, IMXER is set.
1571 C Then Y is loaded from YH and T is set to TN.
1572 C The optional outputs are loaded into the work arrays before returning.
1573 C-----------------------------------------------------------------------
1574 C The maximum number of steps was taken before reaching TOUT. ----------
1581 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1587 C Too much accuracy requested for machine precision. -------------------
1595 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1602 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1611 C KFLAG = -3. Unrecoverable error from PSOL. --------------------------
1618 C Compute IMXER if relevant. -------------------------------------------
1624 BIG = SIZE
1625 IMXER = I
1628 C Set Y vector, T, and optional outputs. -------------------------------
1631 T = TN
1645 RETURN
1646 C-----------------------------------------------------------------------
1647 C Block I.
1648 C The following block handles all error returns due to illegal input
1649 C (ISTATE = -3), as detected before calling the core integrator.
1650 C First the error message routine is called. If the illegal input
1651 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1652 C-----------------------------------------------------------------------
1734 C
1736 RETURN
1737 C
1740 RETURN
1741 C----------------------- End of Subroutine DLSODPK ---------------------
1742 END