| blob | fdfe613c9d2c163d2aa3f13af8b4fa75679f7fce |
1 *DECK DLSODI
9 C-----------------------------------------------------------------------
10 C This is the 18 November 2003 version of
11 C DLSODI: Livermore Solver for Ordinary Differential Equations
12 C (Implicit form).
13 C
14 C This version is in double precision.
15 C
16 C DLSODI solves the initial value problem for linearly implicit
17 C systems of first order ODEs,
18 C A(t,y) * dy/dt = g(t,y) , where A(t,y) is a square matrix,
19 C or, in component form,
20 C ( a * ( dy / dt )) + ... + ( a * ( dy / dt )) =
21 C i,1 1 i,NEQ NEQ
22 C
23 C = g ( t, y , y ,..., y ) ( i = 1,...,NEQ )
24 C i 1 2 NEQ
25 C
26 C If A is singular, this is a differential-algebraic system.
27 C
28 C DLSODI is a variant version of the DLSODE package.
29 C-----------------------------------------------------------------------
30 C Reference:
31 C Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE
32 C Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.),
33 C North-Holland, Amsterdam, 1983, pp. 55-64.
34 C-----------------------------------------------------------------------
35 C Authors: Alan C. Hindmarsh and Jeffrey F. Painter
36 C Center for Applied Scientific Computing, L-561
37 C Lawrence Livermore National Laboratory
38 C Livermore, CA 94551
39 C-----------------------------------------------------------------------
40 C Summary of Usage.
41 C
42 C Communication between the user and the DLSODI package, for normal
43 C situations, is summarized here. This summary describes only a subset
44 C of the full set of options available. See the full description for
45 C details, including optional communication, nonstandard options,
46 C and instructions for special situations. See also the example
47 C problem (with program and output) following this summary.
48 C
49 C A. First, provide a subroutine of the form:
50 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
51 C DOUBLE PRECISION T, Y(*), S(*), R(*)
52 C which computes the residual function
53 C r = g(t,y) - A(t,y) * s ,
54 C as a function of t and the vectors y and s. (s is an internally
55 C generated approximation to dy/dt.) The arrays Y and S are inputs
56 C to the RES routine and should not be altered. The residual
57 C vector is to be stored in the array R. The argument IRES should be
58 C ignored for casual use of DLSODI. (For uses of IRES, see the
59 C paragraph on RES in the full description below.)
60 C
61 C B. Next, decide whether full or banded form is more economical
62 C for the storage of matrices. DLSODI must deal internally with the
63 C matrices A and dr/dy, where r is the residual function defined above.
64 C DLSODI generates a linear combination of these two matrices, and
65 C this is treated in either full or banded form.
66 C The matrix structure is communicated by a method flag MF,
67 C which is 21 or 22 for the full case, and 24 or 25 in the band case.
68 C In the banded case, DLSODI requires two half-bandwidth
69 C parameters ML and MU. These are, respectively, the widths of the
70 C lower and upper parts of the band, excluding the main diagonal.
71 C Thus the band consists of the locations (i,j) with
72 C i-ML .le. j .le. i+MU, and the full bandwidth is ML+MU+1.
73 C Note that the band must accommodate the nonzero elements of
74 C A(t,y), dg/dy, and d(A*s)/dy (s fixed). Alternatively, one
75 C can define a band that encloses only the elements that are relatively
76 C large in magnitude, and gain some economy in storage and possibly
77 C also efficiency, although the appropriate threshhold for
78 C retaining matrix elements is highly problem-dependent.
79 C
80 C C. You must also provide a subroutine of the form:
81 C SUBROUTINE ADDA (NEQ, T, Y, ML, MU, P, NROWP)
82 C DOUBLE PRECISION T, Y(*), P(NROWP,*)
83 C which adds the matrix A = A(t,y) to the contents of the array P.
84 C T and the Y array are input and should not be altered.
85 C In the full matrix case, this routine should add elements of
86 C to P in the usual order. I.e., add A(i,j) to P(i,j). (Ignore the
87 C ML and MU arguments in this case.)
88 C In the band matrix case, this routine should add element A(i,j)
89 C to P(i-j+MU+1,j). I.e., add the diagonal lines of A to the rows of
90 C P from the top down (the top line of A added to the first row of P).
91 C
92 C D. For the sake of efficiency, you are encouraged to supply the
93 C Jacobian matrix dr/dy in closed form, where r = g(t,y) - A(t,y)*s
94 C (s = a fixed vector) as above. If dr/dy is being supplied,
95 C use MF = 21 or 24, and provide a subroutine of the form:
96 C SUBROUTINE JAC (NEQ, T, Y, S, ML, MU, P, NROWP)
97 C DOUBLE PRECISION T, Y(*), S(*), P(NROWP,*)
98 C which computes dr/dy as a function of t, y, and s. Here T, Y, and
99 C S are inputs, and the routine is to load dr/dy into P as follows:
100 C In the full matrix case (MF = 21), load P(i,j) with dr(i)/dy(j),
101 C the partial derivative of r(i) with respect to y(j). (Ignore the
102 C ML and MU arguments in this case.)
103 C In the band matrix case (MF = 24), load P(i-j+mu+1,j) with
104 C dr(i)/dy(j), i.e. load the diagonal lines of dr/dy into the rows of
105 C P from the top down.
106 C In either case, only nonzero elements need be loaded, and the
107 C indexing of P is the same as in the ADDA routine.
108 C Note that if A is independent of y (or this dependence
109 C is weak enough to be ignored) then JAC is to compute dg/dy.
110 C If it is not feasible to provide a JAC routine, use
111 C MF = 22 or 25, and DLSODI will compute an approximate Jacobian
112 C internally by difference quotients.
113 C
114 C E. Next decide whether or not to provide the initial value of the
115 C derivative vector dy/dt. If the initial value of A(t,y) is
116 C nonsingular (and not too ill-conditioned), you may let DLSODI compute
117 C this vector (ISTATE = 0). (DLSODI will solve the system A*s = g for
118 C s, with initial values of A and g.) If A(t,y) is initially
119 C singular, then the system is a differential-algebraic system, and
120 C you must make use of the particular form of the system to compute the
121 C initial values of y and dy/dt. In that case, use ISTATE = 1 and
122 C load the initial value of dy/dt into the array YDOTI.
123 C The input array YDOTI and the initial Y array must be consistent with
124 C the equations A*dy/dt = g. This implies that the initial residual
125 C r = g(t,y) - A(t,y)*YDOTI must be approximately zero.
126 C
127 C F. Write a main program which calls Subroutine DLSODI once for
128 C each point at which answers are desired. This should also provide
129 C for possible use of logical unit 6 for output of error messages
130 C by DLSODI. On the first call to DLSODI, supply arguments as follows:
131 C RES = name of user subroutine for residual function r.
132 C ADDA = name of user subroutine for computing and adding A(t,y).
133 C JAC = name of user subroutine for Jacobian matrix dr/dy
134 C (MF = 21 or 24). If not used, pass a dummy name.
135 C Note: the names for the RES and ADDA routines and (if used) the
136 C JAC routine must be declared External in the calling program.
137 C NEQ = number of scalar equations in the system.
138 C Y = array of initial values, of length NEQ.
139 C YDOTI = array of length NEQ (containing initial dy/dt if ISTATE = 1).
140 C T = the initial value of the independent variable.
141 C TOUT = first point where output is desired (.ne. T).
142 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
143 C RTOL = relative tolerance parameter (scalar).
144 C ATOL = absolute tolerance parameter (scalar or array).
145 C the estimated local error in y(i) will be controlled so as
146 C to be roughly less (in magnitude) than
147 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
148 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
149 C Thus the local error test passes if, in each component,
150 C either the absolute error is less than ATOL (or ATOL(i)),
151 C or the relative error is less than RTOL.
152 C Use RTOL = 0.0 for pure absolute error control, and
153 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
154 C control. Caution: Actual (global) errors may exceed these
155 C local tolerances, so choose them conservatively.
156 C ITASK = 1 for normal computation of output values of y at t = TOUT.
157 C ISTATE = integer flag (input and output). Set ISTATE = 1 if the
158 C initial dy/dt is supplied, and 0 otherwise.
159 C IOPT = 0 to indicate no optional inputs used.
160 C RWORK = real work array of length at least:
161 C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
162 C 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25.
163 C LRW = declared length of RWORK (in user's dimension).
164 C IWORK = integer work array of length at least 20 + NEQ.
165 C If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower
166 C and upper half-bandwidths ML,MU.
167 C LIW = declared length of IWORK (in user's dimension).
168 C MF = method flag. Standard values are:
169 C 21 for a user-supplied full Jacobian.
170 C 22 for an internally generated full Jacobian.
171 C 24 for a user-supplied banded Jacobian.
172 C 25 for an internally generated banded Jacobian.
173 C for other choices of MF, see the paragraph on MF in
174 C the full description below.
175 C Note that the main program must declare arrays Y, YDOTI, RWORK, IWORK,
176 C and possibly ATOL.
177 C
178 C G. The output from the first call (or any call) is:
179 C Y = array of computed values of y(t) vector.
180 C T = corresponding value of independent variable (normally TOUT).
181 C ISTATE = 2 if DLSODI was successful, negative otherwise.
182 C -1 means excess work done on this call (check all inputs).
183 C -2 means excess accuracy requested (tolerances too small).
184 C -3 means illegal input detected (see printed message).
185 C -4 means repeated error test failures (check all inputs).
186 C -5 means repeated convergence failures (perhaps bad Jacobian
187 C supplied or wrong choice of tolerances).
188 C -6 means error weight became zero during problem. (Solution
189 C component i vanished, and ATOL or ATOL(i) = 0.)
190 C -7 cannot occur in casual use.
191 C -8 means DLSODI was unable to compute the initial dy/dt.
192 C In casual use, this means A(t,y) is initially singular.
193 C Supply YDOTI and use ISTATE = 1 on the first call.
194 C
195 C If DLSODI returns ISTATE = -1, -4, or -5, then the output of
196 C DLSODI also includes YDOTI = array containing residual vector
197 C r = g - A * dy/dt evaluated at the current t, y, and dy/dt.
198 C
199 C H. To continue the integration after a successful return, simply
200 C reset TOUT and call DLSODI again. No other parameters need be reset.
201 C
202 C-----------------------------------------------------------------------
203 C Example Problem.
204 C
205 C The following is a simple example problem, with the coding
206 C needed for its solution by DLSODI. The problem is from chemical
207 C kinetics, and consists of the following three equations:
208 C dy1/dt = -.04*y1 + 1.e4*y2*y3
209 C dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
210 C 0. = y1 + y2 + y3 - 1.
211 C on the interval from t = 0.0 to t = 4.e10, with initial conditions
212 C y1 = 1.0, y2 = y3 = 0.
213 C
214 C The following coding solves this problem with DLSODI, using MF = 21
215 C and printing results at t = .4, 4., ..., 4.e10. It uses
216 C ITOL = 2 and ATOL much smaller for y2 than y1 or y3 because
217 C y2 has much smaller values. dy/dt is supplied in YDOTI. We had
218 C obtained the initial value of dy3/dt by differentiating the
219 C third equation and evaluating the first two at t = 0.
220 C At the end of the run, statistical quantities of interest are
221 C printed (see optional outputs in the full description below).
222 C
223 C EXTERNAL RESID, APLUSP, DGBYDY
224 C DOUBLE PRECISION ATOL, RTOL, RWORK, T, TOUT, Y, YDOTI
225 C DIMENSION Y(3), YDOTI(3), ATOL(3), RWORK(58), IWORK(23)
226 C NEQ = 3
227 C Y(1) = 1.
228 C Y(2) = 0.
229 C Y(3) = 0.
230 C YDOTI(1) = -.04
231 C YDOTI(2) = .04
232 C YDOTI(3) = 0.
233 C T = 0.
234 C TOUT = .4
235 C ITOL = 2
236 C RTOL = 1.D-4
237 C ATOL(1) = 1.D-6
238 C ATOL(2) = 1.D-10
239 C ATOL(3) = 1.D-6
240 C ITASK = 1
241 C ISTATE = 1
242 C IOPT = 0
243 C LRW = 58
244 C LIW = 23
245 C MF = 21
246 C DO 40 IOUT = 1,12
247 C CALL DLSODI(RESID, APLUSP, DGBYDY, NEQ, Y, YDOTI, T, TOUT, ITOL,
248 C 1 RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, MF)
249 C WRITE (6,20) T, Y(1), Y(2), Y(3)
250 C 20 FORMAT(' At t =',D12.4,' Y =',3D14.6)
251 C IF (ISTATE .LT. 0 ) GO TO 80
252 C 40 TOUT = TOUT*10.
253 C WRITE (6,60) IWORK(11), IWORK(12), IWORK(13)
254 C 60 FORMAT(/' No. steps =',I4,' No. r-s =',I4,' No. J-s =',I4)
255 C STOP
256 C 80 WRITE (6,90) ISTATE
257 C 90 FORMAT(///' Error halt.. ISTATE =',I3)
258 C STOP
259 C END
260 C
261 C SUBROUTINE RESID(NEQ, T, Y, S, R, IRES)
262 C DOUBLE PRECISION T, Y, S, R
263 C DIMENSION Y(3), S(3), R(3)
264 C R(1) = -.04*Y(1) + 1.D4*Y(2)*Y(3) - S(1)
265 C R(2) = .04*Y(1) - 1.D4*Y(2)*Y(3) - 3.D7*Y(2)*Y(2) - S(2)
266 C R(3) = Y(1) + Y(2) + Y(3) - 1.
267 C RETURN
268 C END
269 C
270 C SUBROUTINE APLUSP(NEQ, T, Y, ML, MU, P, NROWP)
271 C DOUBLE PRECISION T, Y, P
272 C DIMENSION Y(3), P(NROWP,3)
273 C P(1,1) = P(1,1) + 1.
274 C P(2,2) = P(2,2) + 1.
275 C RETURN
276 C END
277 C
278 C SUBROUTINE DGBYDY(NEQ, T, Y, S, ML, MU, P, NROWP)
279 C DOUBLE PRECISION T, Y, S, P
280 C DIMENSION Y(3), S(3), P(NROWP,3)
281 C P(1,1) = -.04
282 C P(1,2) = 1.D4*Y(3)
283 C P(1,3) = 1.D4*Y(2)
284 C P(2,1) = .04
285 C P(2,2) = -1.D4*Y(3) - 6.D7*Y(2)
286 C P(2,3) = -1.D4*Y(2)
287 C P(3,1) = 1.
288 C P(3,2) = 1.
289 C P(3,3) = 1.
290 C RETURN
291 C END
292 C
293 C The output of this program (on a CDC-7600 in single precision)
294 C is as follows:
295 C
296 C At t = 4.0000e-01 Y = 9.851726e-01 3.386406e-05 1.479357e-02
297 C At t = 4.0000e+00 Y = 9.055142e-01 2.240418e-05 9.446344e-02
298 C At t = 4.0000e+01 Y = 7.158050e-01 9.184616e-06 2.841858e-01
299 C At t = 4.0000e+02 Y = 4.504846e-01 3.222434e-06 5.495122e-01
300 C At t = 4.0000e+03 Y = 1.831701e-01 8.940379e-07 8.168290e-01
301 C At t = 4.0000e+04 Y = 3.897016e-02 1.621193e-07 9.610297e-01
302 C At t = 4.0000e+05 Y = 4.935213e-03 1.983756e-08 9.950648e-01
303 C At t = 4.0000e+06 Y = 5.159269e-04 2.064759e-09 9.994841e-01
304 C At t = 4.0000e+07 Y = 5.306413e-05 2.122677e-10 9.999469e-01
305 C At t = 4.0000e+08 Y = 5.494532e-06 2.197826e-11 9.999945e-01
306 C At t = 4.0000e+09 Y = 5.129457e-07 2.051784e-12 9.999995e-01
307 C At t = 4.0000e+10 Y = -7.170472e-08 -2.868188e-13 1.000000e+00
308 C
309 C No. steps = 330 No. r-s = 404 No. J-s = 69
310 C
311 C-----------------------------------------------------------------------
312 C Full Description of User Interface to DLSODI.
313 C
314 C The user interface to DLSODI consists of the following parts.
315 C
316 C 1. The call sequence to Subroutine DLSODI, which is a driver
317 C routine for the solver. This includes descriptions of both
318 C the call sequence arguments and of user-supplied routines.
319 C Following these descriptions is a description of
320 C optional inputs available through the call sequence, and then
321 C a description of optional outputs (in the work arrays).
322 C
323 C 2. Descriptions of other routines in the DLSODI package that may be
324 C (optionally) called by the user. These provide the ability to
325 C alter error message handling, save and restore the internal
326 C Common, and obtain specified derivatives of the solution y(t).
327 C
328 C 3. Descriptions of Common blocks to be declared in overlay
329 C or similar environments, or to be saved when doing an interrupt
330 C of the problem and continued solution later.
331 C
332 C 4. Description of two routines in the DLSODI package, either of
333 C which the user may replace with his/her own version, if desired.
334 C These relate to the measurement of errors.
335 C
336 C-----------------------------------------------------------------------
337 C Part 1. Call Sequence.
338 C
339 C The call sequence parameters used for input only are
340 C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
341 C IOPT, LRW, LIW, MF,
342 C and those used for both input and output are
343 C Y, T, ISTATE, YDOTI.
344 C The work arrays RWORK and IWORK are also used for conditional and
345 C optional inputs and optional outputs. (The term output here refers
346 C to the return from Subroutine DLSODI to the user's calling program.)
347 C
348 C The legality of input parameters will be thoroughly checked on the
349 C initial call for the problem, but not checked thereafter unless a
350 C change in input parameters is flagged by ISTATE = 3 on input.
351 C
352 C The descriptions of the call arguments are as follows.
353 C
354 C RES = the name of the user-supplied subroutine which supplies
355 C the residual vector for the ODE system, defined by
356 C r = g(t,y) - A(t,y) * s
357 C as a function of the scalar t and the vectors
358 C s and y (s approximates dy/dt). This subroutine
359 C is to have the form
360 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
361 C DOUBLE PRECISION T, Y(*), S(*), R(*)
362 C where NEQ, T, Y, S, and IRES are input, and R and
363 C IRES are output. Y, S, and R are arrays of length NEQ.
364 C On input, IRES indicates how DLSODI will use the
365 C returned array R, as follows:
366 C IRES = 1 means that DLSODI needs the full residual,
367 C r = g - A*s, exactly.
368 C IRES = -1 means that DLSODI is using R only to compute
369 C the Jacobian dr/dy by difference quotients.
370 C The RES routine can ignore IRES, or it can omit some terms
371 C if IRES = -1. If A does not depend on y, then RES can
372 C just return R = g when IRES = -1. If g - A*s contains other
373 C additive terms that are independent of y, these can also be
374 C dropped, if done consistently, when IRES = -1.
375 C The subroutine should set the flag IRES if it
376 C encounters a halt condition or illegal input.
377 C Otherwise, it should not reset IRES. On output,
378 C IRES = 1 or -1 represents a normal return, and
379 C DLSODI continues integrating the ODE. Leave IRES
380 C unchanged from its input value.
381 C IRES = 2 tells DLSODI to immediately return control
382 C to the calling program, with ISTATE = 3. This lets
383 C the calling program change parameters of the problem,
384 C if necessary.
385 C IRES = 3 represents an error condition (for example, an
386 C illegal value of y). DLSODI tries to integrate the system
387 C without getting IRES = 3 from RES. If it cannot, DLSODI
388 C returns with ISTATE = -7 or -1.
389 C On an DLSODI return with ISTATE = 3, -1, or -7, the values
390 C of T and Y returned correspond to the last point reached
391 C successfully without getting the flag IRES = 2 or 3.
392 C The flag values IRES = 2 and 3 should not be used to
393 C handle switches or root-stop conditions. This is better
394 C done by calling DLSODI in a one-step mode and checking the
395 C stopping function for a sign change at each step.
396 C If quantities computed in the RES routine are needed
397 C externally to DLSODI, an extra call to RES should be made
398 C for this purpose, for consistent and accurate results.
399 C To get the current dy/dt for the S argument, use DINTDY.
400 C RES must be declared External in the calling
401 C program. See note below for more about RES.
402 C
403 C ADDA = the name of the user-supplied subroutine which adds the
404 C matrix A = A(t,y) to another matrix stored in the same form
405 C as A. The storage form is determined by MITER (see MF).
406 C This subroutine is to have the form
407 C SUBROUTINE ADDA (NEQ, T, Y, ML, MU, P, NROWP)
408 C DOUBLE PRECISION T, Y(*), P(NROWP,*)
409 C where NEQ, T, Y, ML, MU, and NROWP are input and P is
410 C output. Y is an array of length NEQ, and the matrix P is
411 C stored in an NROWP by NEQ array.
412 C In the full matrix case ( MITER = 1 or 2) ADDA should
413 C add A to P(i,j). ML and MU are ignored.
414 C i,j
415 C In the band matrix case ( MITER = 4 or 5) ADDA should
416 C add A to P(i-j+MU+1,j).
417 C i,j
418 C See JAC for details on this band storage form.
419 C ADDA must be declared External in the calling program.
420 C See note below for more information about ADDA.
421 C
422 C JAC = the name of the user-supplied subroutine which supplies the
423 C Jacobian matrix, dr/dy, where r = g - A*s. The form of the
424 C Jacobian matrix is determined by MITER. JAC is required
425 C if MITER = 1 or 4 -- otherwise a dummy name can be
426 C passed. This subroutine is to have the form
427 C SUBROUTINE JAC ( NEQ, T, Y, S, ML, MU, P, NROWP )
428 C DOUBLE PRECISION T, Y(*), S(*), P(NROWP,*)
429 C where NEQ, T, Y, S, ML, MU, and NROWP are input and P
430 C is output. Y and S are arrays of length NEQ, and the
431 C matrix P is stored in an NROWP by NEQ array.
432 C P is to be loaded with partial derivatives (elements
433 C of the Jacobian matrix) on output.
434 C In the full matrix case (MITER = 1), ML and MU
435 C are ignored and the Jacobian is to be loaded into P
436 C by columns-- i.e., dr(i)/dy(j) is loaded into P(i,j).
437 C In the band matrix case (MITER = 4), the elements
438 C within the band are to be loaded into P by columns,
439 C with diagonal lines of dr/dy loaded into the
440 C rows of P. Thus dr(i)/dy(j) is to be loaded
441 C into P(i-j+MU+1,j). The locations in P in the two
442 C triangular areas which correspond to nonexistent matrix
443 C elements can be ignored or loaded arbitrarily, as they
444 C they are overwritten by DLSODI. ML and MU are the
445 C half-bandwidth parameters (see IWORK).
446 C In either case, P is preset to zero by the solver,
447 C so that only the nonzero elements need be loaded by JAC.
448 C Each call to JAC is preceded by a call to RES with the same
449 C arguments NEQ, T, Y, and S. Thus to gain some efficiency,
450 C intermediate quantities shared by both calculations may be
451 C saved in a user Common block by RES and not recomputed by JAC
452 C if desired. Also, JAC may alter the Y array, if desired.
453 C JAC need not provide dr/dy exactly. A crude
454 C approximation (possibly with a smaller bandwidth) will do.
455 C JAC must be declared External in the calling program.
456 C See note below for more about JAC.
457 C
458 C Note on RES, ADDA, and JAC:
459 C These subroutines may access user-defined quantities in
460 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
461 C (dimensioned in the subroutines) and/or Y has length
462 C exceeding NEQ(1). However, these routines should not alter
463 C NEQ(1), Y(1),...,Y(NEQ) or any other input variables.
464 C See the descriptions of NEQ and Y below.
465 C
466 C NEQ = the size of the system (number of first order ordinary
467 C differential equations or scalar algebraic equations).
468 C Used only for input.
469 C NEQ may be decreased, but not increased, during the problem.
470 C If NEQ is decreased (with ISTATE = 3 on input), the
471 C remaining components of Y should be left undisturbed, if
472 C these are to be accessed in RES, ADDA, or JAC.
473 C
474 C Normally, NEQ is a scalar, and it is generally referred to
475 C as a scalar in this user interface description. However,
476 C NEQ may be an array, with NEQ(1) set to the system size.
477 C (The DLSODI package accesses only NEQ(1).) In either case,
478 C this parameter is passed as the NEQ argument in all calls
479 C to RES, ADDA, and JAC. Hence, if it is an array,
480 C locations NEQ(2),... may be used to store other integer data
481 C and pass it to RES, ADDA, or JAC. Each such subroutine
482 C must include NEQ in a Dimension statement in that case.
483 C
484 C Y = a real array for the vector of dependent variables, of
485 C length NEQ or more. Used for both input and output on the
486 C first call (ISTATE = 0 or 1), and only for output on other
487 C calls. On the first call, Y must contain the vector of
488 C initial values. On output, Y contains the computed solution
489 C vector, evaluated at T. If desired, the Y array may be used
490 C for other purposes between calls to the solver.
491 C
492 C This array is passed as the Y argument in all calls to RES,
493 C ADDA, and JAC. Hence its length may exceed NEQ,
494 C and locations Y(NEQ+1),... may be used to store other real
495 C data and pass it to RES, ADDA, or JAC. (The DLSODI
496 C package accesses only Y(1),...,Y(NEQ). )
497 C
498 C YDOTI = a real array for the initial value of the vector
499 C dy/dt and for work space, of dimension at least NEQ.
500 C
501 C On input:
502 C If ISTATE = 0, then DLSODI will compute the initial value
503 C of dy/dt, if A is nonsingular. Thus YDOTI will
504 C serve only as work space and may have any value.
505 C If ISTATE = 1, then YDOTI must contain the initial value
506 C of dy/dt.
507 C If ISTATE = 2 or 3 (continuation calls), then YDOTI
508 C may have any value.
509 C Note: If the initial value of A is singular, then
510 C DLSODI cannot compute the initial value of dy/dt, so
511 C it must be provided in YDOTI, with ISTATE = 1.
512 C
513 C On output, when DLSODI terminates abnormally with ISTATE =
514 C -1, -4, or -5, YDOTI will contain the residual
515 C r = g(t,y) - A(t,y)*(dy/dt). If r is large, t is near
516 C its initial value, and YDOTI is supplied with ISTATE = 1,
517 C then there may have been an incorrect input value of
518 C YDOTI = dy/dt, or the problem (as given to DLSODI)
519 C may not have a solution.
520 C
521 C If desired, the YDOTI array may be used for other
522 C purposes between calls to the solver.
523 C
524 C T = the independent variable. On input, T is used only on the
525 C first call, as the initial point of the integration.
526 C On output, after each call, T is the value at which a
527 C computed solution Y is evaluated (usually the same as TOUT).
528 C on an error return, T is the farthest point reached.
529 C
530 C TOUT = the next value of t at which a computed solution is desired.
531 C Used only for input.
532 C
533 C When starting the problem (ISTATE = 0 or 1), TOUT may be
534 C equal to T for one call, then should .ne. T for the next
535 C call. For the initial T, an input value of TOUT .ne. T is
536 C used in order to determine the direction of the integration
537 C (i.e. the algebraic sign of the step sizes) and the rough
538 C scale of the problem. Integration in either direction
539 C (forward or backward in t) is permitted.
540 C
541 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
542 C the first call (i.e. the first call with TOUT .ne. T).
543 C Otherwise, TOUT is required on every call.
544 C
545 C If ITASK = 1, 3, or 4, the values of TOUT need not be
546 C monotone, but a value of TOUT which backs up is limited
547 C to the current internal T interval, whose endpoints are
548 C TCUR - HU and TCUR (see optional outputs, below, for
549 C TCUR and HU).
550 C
551 C ITOL = an indicator for the type of error control. See
552 C description below under ATOL. Used only for input.
553 C
554 C RTOL = a relative error tolerance parameter, either a scalar or
555 C an array of length NEQ. See description below under ATOL.
556 C Input only.
557 C
558 C ATOL = an absolute error tolerance parameter, either a scalar or
559 C an array of length NEQ. Input only.
560 C
561 C The input parameters ITOL, RTOL, and ATOL determine
562 C the error control performed by the solver. The solver will
563 C control the vector E = (E(i)) of estimated local errors
564 C in y, according to an inequality of the form
565 C RMS-norm of ( E(i)/EWT(i) ) .le. 1,
566 C where EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
567 C and the RMS-norm (root-mean-square norm) here is
568 C RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
569 C is a vector of weights which must always be positive, and
570 C the values of RTOL and ATOL should all be non-negative.
571 C The following table gives the types (scalar/array) of
572 C RTOL and ATOL, and the corresponding form of EWT(i).
573 C
574 C ITOL RTOL ATOL EWT(i)
575 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
576 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
577 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
578 C 4 array scalar RTOL(i)*ABS(Y(i)) + ATOL(i)
579 C
580 C When either of these parameters is a scalar, it need not
581 C be dimensioned in the user's calling program.
582 C
583 C If none of the above choices (with ITOL, RTOL, and ATOL
584 C fixed throughout the problem) is suitable, more general
585 C error controls can be obtained by substituting
586 C user-supplied routines for the setting of EWT and/or for
587 C the norm calculation. See Part 4 below.
588 C
589 C If global errors are to be estimated by making a repeated
590 C run on the same problem with smaller tolerances, then all
591 C components of RTOL and ATOL (i.e. of EWT) should be scaled
592 C down uniformly.
593 C
594 C ITASK = an index specifying the task to be performed.
595 C Input only. ITASK has the following values and meanings.
596 C 1 means normal computation of output values of y(t) at
597 C t = TOUT (by overshooting and interpolating).
598 C 2 means take one step only and return.
599 C 3 means stop at the first internal mesh point at or
600 C beyond t = TOUT and return.
601 C 4 means normal computation of output values of y(t) at
602 C t = TOUT but without overshooting t = TCRIT.
603 C TCRIT must be input as RWORK(1). TCRIT may be equal to
604 C or beyond TOUT, but not behind it in the direction of
605 C integration. This option is useful if the problem
606 C has a singularity at or beyond t = TCRIT.
607 C 5 means take one step, without passing TCRIT, and return.
608 C TCRIT must be input as RWORK(1).
609 C
610 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
611 C (within roundoff), it will return T = TCRIT (exactly) to
612 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
613 C in which case answers at t = TOUT are returned first).
614 C
615 C ISTATE = an index used for input and output to specify the
616 C state of the calculation.
617 C
618 C On input, the values of ISTATE are as follows.
619 C 0 means this is the first call for the problem, and
620 C DLSODI is to compute the initial value of dy/dt
621 C (while doing other initializations). See note below.
622 C 1 means this is the first call for the problem, and
623 C the initial value of dy/dt has been supplied in
624 C YDOTI (DLSODI will do other initializations). See note
625 C below.
626 C 2 means this is not the first call, and the calculation
627 C is to continue normally, with no change in any input
628 C parameters except possibly TOUT and ITASK.
629 C (If ITOL, RTOL, and/or ATOL are changed between calls
630 C with ISTATE = 2, the new values will be used but not
631 C tested for legality.)
632 C 3 means this is not the first call, and the
633 C calculation is to continue normally, but with
634 C a change in input parameters other than
635 C TOUT and ITASK. Changes are allowed in
636 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU,
637 C and any of the optional inputs except H0.
638 C (See IWORK description for ML and MU.)
639 C Note: A preliminary call with TOUT = T is not counted
640 C as a first call here, as no initialization or checking of
641 C input is done. (Such a call is sometimes useful for the
642 C purpose of outputting the initial conditions.)
643 C Thus the first call for which TOUT .ne. T requires
644 C ISTATE = 0 or 1 on input.
645 C
646 C On output, ISTATE has the following values and meanings.
647 C 0 or 1 means nothing was done; TOUT = t and
648 C ISTATE = 0 or 1 on input.
649 C 2 means that the integration was performed successfully.
650 C 3 means that the user-supplied Subroutine RES signalled
651 C DLSODI to halt the integration and return (IRES = 2).
652 C Integration as far as T was achieved with no occurrence
653 C of IRES = 2, but this flag was set on attempting the
654 C next step.
655 C -1 means an excessive amount of work (more than MXSTEP
656 C steps) was done on this call, before completing the
657 C requested task, but the integration was otherwise
658 C successful as far as T. (MXSTEP is an optional input
659 C and is normally 500.) To continue, the user may
660 C simply reset ISTATE to a value .gt. 1 and call again
661 C (the excess work step counter will be reset to 0).
662 C In addition, the user may increase MXSTEP to avoid
663 C this error return (see below on optional inputs).
664 C -2 means too much accuracy was requested for the precision
665 C of the machine being used. This was detected before
666 C completing the requested task, but the integration
667 C was successful as far as T. To continue, the tolerance
668 C parameters must be reset, and ISTATE must be set
669 C to 3. The optional output TOLSF may be used for this
670 C purpose. (Note: If this condition is detected before
671 C taking any steps, then an illegal input return
672 C (ISTATE = -3) occurs instead.)
673 C -3 means illegal input was detected, before taking any
674 C integration steps. See written message for details.
675 C Note: If the solver detects an infinite loop of calls
676 C to the solver with illegal input, it will cause
677 C the run to stop.
678 C -4 means there were repeated error test failures on
679 C one attempted step, before completing the requested
680 C task, but the integration was successful as far as T.
681 C The problem may have a singularity, or the input
682 C may be inappropriate.
683 C -5 means there were repeated convergence test failures on
684 C one attempted step, before completing the requested
685 C task, but the integration was successful as far as T.
686 C This may be caused by an inaccurate Jacobian matrix.
687 C -6 means EWT(i) became zero for some i during the
688 C integration. pure relative error control (ATOL(i)=0.0)
689 C was requested on a variable which has now vanished.
690 C the integration was successful as far as T.
691 C -7 means that the user-supplied Subroutine RES set
692 C its error flag (IRES = 3) despite repeated tries by
693 C DLSODI to avoid that condition.
694 C -8 means that ISTATE was 0 on input but DLSODI was unable
695 C to compute the initial value of dy/dt. See the
696 C printed message for details.
697 C
698 C Note: Since the normal output value of ISTATE is 2,
699 C it does not need to be reset for normal continuation.
700 C Similarly, ISTATE (= 3) need not be reset if RES told
701 C DLSODI to return because the calling program must change
702 C the parameters of the problem.
703 C Also, since a negative input value of ISTATE will be
704 C regarded as illegal, a negative output value requires the
705 C user to change it, and possibly other inputs, before
706 C calling the solver again.
707 C
708 C IOPT = an integer flag to specify whether or not any optional
709 C inputs are being used on this call. Input only.
710 C The optional inputs are listed separately below.
711 C IOPT = 0 means no optional inputs are being used.
712 C Default values will be used in all cases.
713 C IOPT = 1 means one or more optional inputs are being used.
714 C
715 C RWORK = a real working array (double precision).
716 C The length of RWORK must be at least
717 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LENWM where
718 C NYH = the initial value of NEQ,
719 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
720 C smaller value is given as an optional input),
721 C LENWM = NEQ**2 + 2 if MITER is 1 or 2, and
722 C LENWM = (2*ML+MU+1)*NEQ + 2 if MITER is 4 or 5.
723 C (See MF description for the definition of METH and MITER.)
724 C Thus if MAXORD has its default value and NEQ is constant,
725 C this length is
726 C 22 + 16*NEQ + NEQ**2 for MF = 11 or 12,
727 C 22 + 17*NEQ + (2*ML+MU)*NEQ for MF = 14 or 15,
728 C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
729 C 22 + 10*NEQ + (2*ML+MU)*NEQ for MF = 24 or 25.
730 C The first 20 words of RWORK are reserved for conditional
731 C and optional inputs and optional outputs.
732 C
733 C The following word in RWORK is a conditional input:
734 C RWORK(1) = TCRIT = critical value of t which the solver
735 C is not to overshoot. Required if ITASK is
736 C 4 or 5, and ignored otherwise. (See ITASK.)
737 C
738 C LRW = the length of the array RWORK, as declared by the user.
739 C (This will be checked by the solver.)
740 C
741 C IWORK = an integer work array. The length of IWORK must be at least
742 C 20 + NEQ . The first few words of IWORK are used for
743 C conditional and optional inputs and optional outputs.
744 C
745 C The following 2 words in IWORK are conditional inputs:
746 C IWORK(1) = ML These are the lower and upper
747 C IWORK(2) = MU half-bandwidths, respectively, of the
748 C matrices in the problem-- the Jacobian dr/dy
749 C and the left-hand side matrix A. These
750 C half-bandwidths exclude the main diagonal,
751 C so the total bandwidth is ML + MU + 1 .
752 C The band is defined by the matrix locations
753 C (i,j) with i-ML .le. j .le. i+MU. ML and MU
754 C must satisfy 0 .le. ML,MU .le. NEQ-1.
755 C These are required if MITER is 4 or 5, and
756 C ignored otherwise.
757 C ML and MU may in fact be the band parameters
758 C for matrices to which dr/dy and A are only
759 C approximately equal.
760 C
761 C LIW = the length of the array IWORK, as declared by the user.
762 C (This will be checked by the solver.)
763 C
764 C Note: The work arrays must not be altered between calls to DLSODI
765 C for the same problem, except possibly for the conditional and
766 C optional inputs, and except for the last 3*NEQ words of RWORK.
767 C The latter space is used for internal scratch space, and so is
768 C available for use by the user outside DLSODI between calls, if
769 C desired (but not for use by RES, ADDA, or JAC).
770 C
771 C MF = the method flag. Used only for input. The legal values of
772 C MF are 11, 12, 14, 15, 21, 22, 24, and 25.
773 C MF has decimal digits METH and MITER: MF = 10*METH + MITER.
774 C METH indicates the basic linear multistep method:
775 C METH = 1 means the implicit Adams method.
776 C METH = 2 means the method based on Backward
777 C Differentiation Formulas (BDFs).
778 C The BDF method is strongly preferred for stiff
779 C problems, while the Adams method is preferred when
780 C the problem is not stiff. If the matrix A(t,y) is
781 C nonsingular, stiffness here can be taken to mean that of
782 C the explicit ODE system dy/dt = A-inverse * g. If A is
783 C singular, the concept of stiffness is not well defined.
784 C If you do not know whether the problem is stiff, we
785 C recommend using METH = 2. If it is stiff, the advantage
786 C of METH = 2 over METH = 1 will be great, while if it is
787 C not stiff, the advantage of METH = 1 will be slight.
788 C If maximum efficiency is important, some experimentation
789 C with METH may be necessary.
790 C MITER indicates the corrector iteration method:
791 C MITER = 1 means chord iteration with a user-supplied
792 C full (NEQ by NEQ) Jacobian.
793 C MITER = 2 means chord iteration with an internally
794 C generated (difference quotient) full Jacobian.
795 C This uses NEQ+1 extra calls to RES per dr/dy
796 C evaluation.
797 C MITER = 4 means chord iteration with a user-supplied
798 C banded Jacobian.
799 C MITER = 5 means chord iteration with an internally
800 C generated banded Jacobian (using ML+MU+2
801 C extra calls to RES per dr/dy evaluation).
802 C If MITER = 1 or 4, the user must supply a Subroutine JAC
803 C (the name is arbitrary) as described above under JAC.
804 C For other values of MITER, a dummy argument can be used.
805 C-----------------------------------------------------------------------
806 C Optional Inputs.
807 C
808 C The following is a list of the optional inputs provided for in the
809 C call sequence. (See also Part 2.) For each such input variable,
810 C this table lists its name as used in this documentation, its
811 C location in the call sequence, its meaning, and the default value.
812 C the use of any of these inputs requires IOPT = 1, and in that
813 C case all of these inputs are examined. A value of zero for any
814 C of these optional inputs will cause the default value to be used.
815 C Thus to use a subset of the optional inputs, simply preload
816 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
817 C then set those of interest to nonzero values.
818 C
819 C Name Location Meaning and Default Value
820 C
821 C H0 RWORK(5) the step size to be attempted on the first step.
822 C The default value is determined by the solver.
823 C
824 C HMAX RWORK(6) the maximum absolute step size allowed.
825 C The default value is infinite.
826 C
827 C HMIN RWORK(7) the minimum absolute step size allowed.
828 C The default value is 0. (This lower bound is not
829 C enforced on the final step before reaching TCRIT
830 C when ITASK = 4 or 5.)
831 C
832 C MAXORD IWORK(5) the maximum order to be allowed. The default
833 C value is 12 if METH = 1, and 5 if METH = 2.
834 C If MAXORD exceeds the default value, it will
835 C be reduced to the default value.
836 C If MAXORD is changed during the problem, it may
837 C cause the current order to be reduced.
838 C
839 C MXSTEP IWORK(6) maximum number of (internally defined) steps
840 C allowed during one call to the solver.
841 C The default value is 500.
842 C
843 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
844 C warning that T + H = T on a step (H = step size).
845 C This must be positive to result in a non-default
846 C value. The default value is 10.
847 C-----------------------------------------------------------------------
848 C Optional Outputs.
849 C
850 C As optional additional output from DLSODI, the variables listed
851 C below are quantities related to the performance of DLSODI
852 C which are available to the user. These are communicated by way of
853 C the work arrays, but also have internal mnemonic names as shown.
854 C Except where stated otherwise, all of these outputs are defined
855 C on any successful return from DLSODI, and on any return with
856 C ISTATE = -1, -2, -4, -5, -6, or -7. On a return with -3 (illegal
857 C input) or -8, they will be unchanged from their existing values
858 C (if any), except possibly for TOLSF, LENRW, and LENIW.
859 C On any error return, outputs relevant to the error will be defined,
860 C as noted below.
861 C
862 C Name Location Meaning
863 C
864 C HU RWORK(11) the step size in t last used (successfully).
865 C
866 C HCUR RWORK(12) the step size to be attempted on the next step.
867 C
868 C TCUR RWORK(13) the current value of the independent variable
869 C which the solver has actually reached, i.e. the
870 C current internal mesh point in t. On output, TCUR
871 C will always be at least as far as the argument
872 C T, but may be farther (if interpolation was done).
873 C
874 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
875 C computed when a request for too much accuracy was
876 C detected (ISTATE = -3 if detected at the start of
877 C the problem, ISTATE = -2 otherwise). If ITOL is
878 C left unaltered but RTOL and ATOL are uniformly
879 C scaled up by a factor of TOLSF for the next call,
880 C then the solver is deemed likely to succeed.
881 C (The user may also ignore TOLSF and alter the
882 C tolerance parameters in any other way appropriate.)
883 C
884 C NST IWORK(11) the number of steps taken for the problem so far.
885 C
886 C NRE IWORK(12) the number of residual evaluations (RES calls)
887 C for the problem so far.
888 C
889 C NJE IWORK(13) the number of Jacobian evaluations (each involving
890 C an evaluation of A and dr/dy) for the problem so
891 C far. This equals the number of calls to ADDA and
892 C (if MITER = 1 or 4) JAC, and the number of matrix
893 C LU decompositions.
894 C
895 C NQU IWORK(14) the method order last used (successfully).
896 C
897 C NQCUR IWORK(15) the order to be attempted on the next step.
898 C
899 C IMXER IWORK(16) the index of the component of largest magnitude in
900 C the weighted local error vector ( E(i)/EWT(i) ),
901 C on an error return with ISTATE = -4 or -5.
902 C
903 C LENRW IWORK(17) the length of RWORK actually required.
904 C This is defined on normal returns and on an illegal
905 C input return for insufficient storage.
906 C
907 C LENIW IWORK(18) the length of IWORK actually required.
908 C This is defined on normal returns and on an illegal
909 C input return for insufficient storage.
910 C
911 C
912 C The following two arrays are segments of the RWORK array which
913 C may also be of interest to the user as optional outputs.
914 C For each array, the table below gives its internal name,
915 C its base address in RWORK, and its description.
916 C
917 C Name Base Address Description
918 C
919 C YH 21 the Nordsieck history array, of size NYH by
920 C (NQCUR + 1), where NYH is the initial value
921 C of NEQ. For j = 0,1,...,NQCUR, column j+1
922 C of YH contains HCUR**j/factorial(j) times
923 C the j-th derivative of the interpolating
924 C polynomial currently representing the solution,
925 C evaluated at t = TCUR.
926 C
927 C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
928 C corrections on each step, scaled on output to
929 C represent the estimated local error in y on the
930 C last step. This is the vector E in the descrip-
931 C tion of the error control. It is defined only
932 C on a return from DLSODI with ISTATE = 2.
933 C
934 C-----------------------------------------------------------------------
935 C Part 2. Other Routines Callable.
936 C
937 C The following are optional calls which the user may make to
938 C gain additional capabilities in conjunction with DLSODI.
939 C (The routines XSETUN and XSETF are designed to conform to the
940 C SLATEC error handling package.)
941 C
942 C Form of Call Function
943 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
944 C output of messages from DLSODI, if
945 C the default is not desired.
946 C The default value of LUN is 6.
947 C
948 C CALL XSETF(MFLAG) Set a flag to control the printing of
949 C messages by DLSODI.
950 C MFLAG = 0 means do not print. (Danger:
951 C This risks losing valuable information.)
952 C MFLAG = 1 means print (the default).
953 C
954 C Either of the above calls may be made at
955 C any time and will take effect immediately.
956 C
957 C CALL DSRCOM(RSAV,ISAV,JOB) saves and restores the contents of
958 C the internal Common blocks used by
959 C DLSODI (see Part 3 below).
960 C RSAV must be a real array of length 218
961 C or more, and ISAV must be an integer
962 C array of length 37 or more.
963 C JOB=1 means save Common into RSAV/ISAV.
964 C JOB=2 means restore Common from RSAV/ISAV.
965 C DSRCOM is useful if one is
966 C interrupting a run and restarting
967 C later, or alternating between two or
968 C more problems solved with DLSODI.
969 C
970 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
971 C (see below) orders, at a specified point t, if
972 C desired. It may be called only after
973 C a successful return from DLSODI.
974 C
975 C The detailed instructions for using DINTDY are as follows.
976 C The form of the call is:
977 C
978 C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
979 C
980 C The input parameters are:
981 C
982 C T = value of independent variable where answers are desired
983 C (normally the same as the T last returned by DLSODI).
984 C For valid results, T must lie between TCUR - HU and TCUR.
985 C (See optional outputs for TCUR and HU.)
986 C K = integer order of the derivative desired. K must satisfy
987 C 0 .le. K .le. NQCUR, where NQCUR is the current order
988 C (see optional outputs). The capability corresponding
989 C to K = 0, i.e. computing y(T), is already provided
990 C by DLSODI directly. Since NQCUR .ge. 1, the first
991 C derivative dy/dt is always available with DINTDY.
992 C RWORK(21) = the base address of the history array YH.
993 C NYH = column length of YH, equal to the initial value of NEQ.
994 C
995 C The output parameters are:
996 C
997 C DKY = a real array of length NEQ containing the computed value
998 C of the K-th derivative of y(t).
999 C IFLAG = integer flag, returned as 0 if K and T were legal,
1000 C -1 if K was illegal, and -2 if T was illegal.
1001 C On an error return, a message is also written.
1002 C-----------------------------------------------------------------------
1003 C Part 3. Common Blocks.
1004 C
1005 C If DLSODI is to be used in an overlay situation, the user
1006 C must declare, in the primary overlay, the variables in:
1007 C (1) the call sequence to DLSODI, and
1008 C (2) the internal Common block
1009 C /DLS001/ of length 255 (218 double precision words
1010 C followed by 37 integer words),
1011 C
1012 C If DLSODI is used on a system in which the contents of internal
1013 C Common blocks are not preserved between calls, the user should
1014 C declare the above Common block in the calling program to insure
1015 C that their contents are preserved.
1016 C
1017 C If the solution of a given problem by DLSODI is to be interrupted
1018 C and then later continued, such as when restarting an interrupted run
1019 C or alternating between two or more problems, the user should save,
1020 C following the return from the last DLSODI call prior to the
1021 C interruption, the contents of the call sequence variables and the
1022 C internal Common blocks, and later restore these values before the
1023 C next DLSODI call for that problem. To save and restore the Common
1024 C blocks, use Subroutine DSRCOM (see Part 2 above).
1025 C
1026 C-----------------------------------------------------------------------
1027 C Part 4. Optionally Replaceable Solver Routines.
1028 C
1029 C Below are descriptions of two routines in the DLSODI package which
1030 C relate to the measurement of errors. Either routine can be
1031 C replaced by a user-supplied version, if desired. However, since such
1032 C a replacement may have a major impact on performance, it should be
1033 C done only when absolutely necessary, and only with great caution.
1034 C (Note: The means by which the package version of a routine is
1035 C superseded by the user's version may be system-dependent.)
1036 C
1037 C (a) DEWSET.
1038 C The following subroutine is called just before each internal
1039 C integration step, and sets the array of error weights, EWT, as
1040 C described under ITOL/RTOL/ATOL above:
1041 C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1042 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODI call sequence,
1043 C YCUR contains the current dependent variable vector, and
1044 C EWT is the array of weights set by DEWSET.
1045 C
1046 C If the user supplies this subroutine, it must return in EWT(i)
1047 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1048 C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
1049 C routine (see below), and also used by DLSODI in the computation
1050 C of the optional output IMXER, the diagonal Jacobian approximation,
1051 C and the increments for difference quotient Jacobians.
1052 C
1053 C In the user-supplied version of DEWSET, it may be desirable to use
1054 C the current values of derivatives of y. Derivatives up to order NQ
1055 C are available from the history array YH, described above under
1056 C optional outputs. In DEWSET, YH is identical to the YCUR array,
1057 C extended to NQ + 1 columns with a column length of NYH and scale
1058 C factors of H**j/factorial(j). On the first call for the problem,
1059 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1060 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
1061 C can be obtained by including in DEWSET the statements:
1062 C DOUBLE PRECISION RLS
1063 C COMMON /DLS001/ RLS(218),ILS(37)
1064 C NQ = ILS(33)
1065 C NST = ILS(34)
1066 C H = RLS(212)
1067 C Thus, for example, the current value of dy/dt can be obtained as
1068 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
1069 C unnecessary when NST = 0).
1070 C
1071 C (b) DVNORM.
1072 C The following is a real function routine which computes the weighted
1073 C root-mean-square norm of a vector v:
1074 C D = DVNORM (N, V, W)
1075 C where:
1076 C N = the length of the vector,
1077 C V = real array of length N containing the vector,
1078 C W = real array of length N containing weights,
1079 C D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).
1080 C DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
1081 C EWT is as set by Subroutine DEWSET.
1082 C
1083 C If the user supplies this function, it should return a non-negative
1084 C value of DVNORM suitable for use in the error control in DLSODI.
1085 C None of the arguments should be altered by DVNORM.
1086 C For example, a user-supplied DVNORM routine might:
1087 C -substitute a max-norm of (V(i)*W(i)) for the RMS-norm, or
1088 C -ignore some components of V in the norm, with the effect of
1089 C suppressing the error control on those components of y.
1090 C-----------------------------------------------------------------------
1091 C
1092 C***REVISION HISTORY (YYYYMMDD)
1093 C 19800424 DATE WRITTEN
1094 C 19800519 Corrected access of YH on forced order reduction;
1095 C numerous corrections to prologues and other comments.
1096 C 19800617 In main driver, added loading of SQRT(UROUND) in RWORK;
1097 C minor corrections to main prologue.
1098 C 19800903 Corrected ISTATE logic; minor changes in prologue.
1099 C 19800923 Added zero initialization of HU and NQU.
1100 C 19801028 Reorganized RES calls in AINVG, STODI, and PREPJI;
1101 C in LSODI, corrected NRE increment and reset LDY0 at 580;
1102 C numerous corrections to main prologue.
1103 C 19801218 Revised XERRWD routine; minor corrections to main prologue.
1104 C 19810330 Added Common block /LSI001/; use LSODE's INTDY and SOLSY;
1105 C minor corrections to XERRWD and error message at 604;
1106 C minor corrections to declarations; corrections to prologues.
1107 C 19810818 Numerous revisions: replaced EWT by 1/EWT; used flags
1108 C JCUR, ICF, IERPJ, IERSL between STODI and subordinates;
1109 C added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF;
1110 C reorganized returns from STODI; reorganized type decls.;
1111 C fixed message length in XERRWD; changed default LUNIT to 6;
1112 C changed Common lengths; changed comments throughout.
1113 C 19820906 Corrected use of ABS(H) in STODI; minor comment fixes.
1114 C 19830510 Numerous revisions: revised diff. quotient increment;
1115 C eliminated block /LSI001/, using IERPJ flag;
1116 C revised STODI logic after PJAC return;
1117 C revised tuning of H change and step attempts in STODI;
1118 C corrections to main prologue and internal comments.
1119 C 19870330 Major update: corrected comments throughout;
1120 C removed TRET from Common; rewrote EWSET with 4 loops;
1121 C fixed t test in INTDY; added Cray directives in STODI;
1122 C in STODI, fixed DELP init. and logic around PJAC call;
1123 C combined routines to save/restore Common;
1124 C passed LEVEL = 0 in error message calls (except run abort).
1125 C 20010425 Major update: convert source lines to upper case;
1126 C added *DECK lines; changed from 1 to * in dummy dimensions;
1127 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
1128 C renamed routines for uniqueness across single/double prec.;
1129 C converted intrinsic names to generic form;
1130 C removed ILLIN and NTREP (data loaded) from Common;
1131 C removed all 'own' variables from Common;
1132 C changed error messages to quoted strings;
1133 C replaced XERRWV/XERRWD with 1993 revised version;
1134 C converted prologues, comments, error messages to mixed case;
1135 C converted arithmetic IF statements to logical IF statements;
1136 C numerous corrections to prologues and internal comments.
1137 C 20010507 Converted single precision source to double precision.
1138 C 20020502 Corrected declarations in descriptions of user routines.
1139 C 20031105 Restored 'own' variables to Common block, to enable
1140 C interrupt/restart feature.
1141 C 20031112 Added SAVE statements for data-loaded constants.
1142 C 20031117 Changed internal names NRE, LSAVR to NFE, LSAVF resp.
1143 C
1144 C-----------------------------------------------------------------------
1145 C Other routines in the DLSODI package.
1146 C
1147 C In addition to Subroutine DLSODI, the DLSODI package includes the
1148 C following subroutines and function routines:
1149 C DAINVG computes the initial value of the vector
1150 C dy/dt = A-inverse * g
1151 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1152 C DSTODI is the core integrator, which does one step of the
1153 C integration and the associated error control.
1154 C DCFODE sets all method coefficients and test constants.
1155 C DPREPJI computes and preprocesses the Jacobian matrix
1156 C and the Newton iteration matrix P.
1157 C DSOLSY manages solution of linear system in chord iteration.
1158 C DEWSET sets the error weight vector EWT before each step.
1159 C DVNORM computes the weighted RMS-norm of a vector.
1160 C DSRCOM is a user-callable routine to save and restore
1161 C the contents of the internal Common blocks.
1162 C DGEFA and DGESL are routines from LINPACK for solving full
1163 C systems of linear algebraic equations.
1164 C DGBFA and DGBSL are routines from LINPACK for solving banded
1165 C linear systems.
1166 C DUMACH computes the unit roundoff in a machine-independent manner.
1167 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
1168 C error messages and warnings. XERRWD is machine-dependent.
1169 C Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1170 C All the others are subroutines.
1171 C
1172 C-----------------------------------------------------------------------
1186 LOGICAL IHIT
1187 CHARACTER*60 MSG
1189 C-----------------------------------------------------------------------
1190 C The following internal Common block contains
1191 C (a) variables which are local to any subroutine but whose values must
1192 C be preserved between calls to the routine ("own" variables), and
1193 C (b) variables which are communicated between subroutines.
1194 C The block DLS001 is declared in subroutines DLSODI, DINTDY, DSTODI,
1195 C DPREPJI, and DSOLSY.
1196 C Groups of variables are replaced by dummy arrays in the Common
1197 C declarations in routines where those variables are not used.
1198 C-----------------------------------------------------------------------
1205 C
1207 C-----------------------------------------------------------------------
1208 C Block A.
1209 C This code block is executed on every call.
1210 C It tests ISTATE and ITASK for legality and branches appropriately.
1211 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1212 C not yet been done, an error return occurs.
1213 C If ISTATE = 0 or 1 and TOUT = T, return immediately.
1214 C-----------------------------------------------------------------------
1223 C-----------------------------------------------------------------------
1224 C Block B.
1225 C The next code block is executed for the initial call (ISTATE = 0 or 1)
1226 C or for a continuation call with parameter changes (ISTATE = 3).
1227 C It contains checking of all inputs and various initializations.
1228 C
1229 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1230 C MF, ML, and MU.
1231 C-----------------------------------------------------------------------
1249 C Next process and check the optional inputs. --------------------------
1252 MXSTEP = MXSTP0
1253 MXHNIL = MXHNL0
1277 C-----------------------------------------------------------------------
1278 C Set work array pointers and check lengths LRW and LIW.
1279 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1280 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1281 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVR, ACOR.
1282 C-----------------------------------------------------------------------
1298 C Check RTOL and ATOL for legality. ------------------------------------
1308 C If ISTATE = 3, set flag to signal parameter changes to DSTODI. -------
1311 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into YDOTI.---------
1314 C Reload WM(1) = RWORK(lWM), since lWM may have changed. ---------------
1317 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1324 C-----------------------------------------------------------------------
1325 C Block C.
1326 C The next block is for the initial call only (ISTATE = 0 or 1).
1327 C It contains all remaining initializations, the call to DAINVG
1328 C (if ISTATE = 1), and the calculation of the initial step size.
1329 C The error weights in EWT are inverted after being loaded.
1330 C-----------------------------------------------------------------------
1332 TN = T
1351 C Compute initial dy/dt, if necessary, and load it and initial Y into YH
1355 C DLSODI must compute initial dy/dt (LYD0 points to YH(*,2)). ----------
1364 C Initial dy/dt was supplied. Load into YH (LYD0 points to YH(*,2).). -
1368 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1376 C-----------------------------------------------------------------------
1377 C The coding below computes the step size, H0, to be attempted on the
1378 C first step, unless the user has supplied a value for this.
1379 C First check that TOUT - T differs significantly from zero.
1380 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1381 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1382 C so as to be between 100*UROUND and 1.0E-3.
1383 C Then the computed value H0 is given by..
1384 C NEQ
1385 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * Sum ( YDOT(i)/ywt(i) )**2 )
1386 C 1
1387 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1388 C YDOT(i) = i-th component of initial value of dy/dt,
1389 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1390 C The sign of H0 is inferred from the initial values of TOUT and T.
1391 C-----------------------------------------------------------------------
1414 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1417 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1418 H = H0
1422 C-----------------------------------------------------------------------
1423 C Block D.
1424 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1425 C and is to check stop conditions before taking a step.
1426 C-----------------------------------------------------------------------
1432 T = TOUT
1444 T = TOUT
1455 C-----------------------------------------------------------------------
1456 C Block E.
1457 C The next block is normally executed for all calls and contains
1458 C the call to the one-step core integrator DSTODI.
1459 C
1460 C This is a looping point for the integration steps.
1461 C
1462 C First check for too many steps being taken, update EWT (if not at
1463 C start of problem), check for too much accuracy being requested, and
1464 C check for H below the roundoff level in T.
1465 C-----------------------------------------------------------------------
1492 C-----------------------------------------------------------------------
1493 C CALL DSTODI(NEQ,Y,YH,NYH,YH1,EWT,SAVF,SAVR,ACOR,WM,IWM,RES,
1494 C ADDA,JAC,DPREPJI,DSOLSY)
1495 C Note: SAVF in DSTODI occupies the same space as YDOTI in DLSODI.
1496 C-----------------------------------------------------------------------
1502 C
1503 C KGO = 1:success; 2:error test failure; 3:convergence failure;
1504 C 4:RES ordered return. 5:RES returned error.
1505 C-----------------------------------------------------------------------
1506 C Block F.
1507 C The following block handles the case of a successful return from the
1508 C core integrator (KFLAG = 0). Test for stop conditions.
1509 C-----------------------------------------------------------------------
1512 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1515 T = TOUT
1517 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1520 C ITASK = 4. see if TOUT or TCRIT was reached. adjust h if necessary.
1523 T = TOUT
1533 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1536 C-----------------------------------------------------------------------
1537 C Block G.
1538 C The following block handles all successful returns from DLSODI.
1539 C if ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1540 C ISTATE is set to 2, and the optional outputs are loaded into the
1541 C work arrays before returning.
1542 C-----------------------------------------------------------------------
1545 T = TN
1558 RETURN
1559 C-----------------------------------------------------------------------
1560 C Block H.
1561 C The following block handles all unsuccessful returns other than
1562 C those for illegal input. First the error message routine is called.
1563 C If there was an error test or convergence test failure, IMXER is set.
1564 C Then Y is loaded from YH and T is set to TN.
1565 C The optional outputs are loaded into the work arrays before returning.
1566 C-----------------------------------------------------------------------
1567 C The maximum number of steps was taken before reaching TOUT. ----------
1574 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1580 C Too much accuracy requested for machine precision. -------------------
1588 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1595 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1604 C IRES = 3 returned by RES, despite retries by DSTODI. -----------------
1611 C DAINVG failed because matrix A was singular. -------------------------
1618 RETURN
1619 C DAINVG failed because RES set IRES to 2 or 3. ------------------------
1627 RETURN
1628 C Compute IMXER if relevant. -------------------------------------------
1634 BIG = SIZE
1635 IMXER = I
1638 C Compute residual if relevant. ----------------------------------------
1652 C Set Y vector, T, and optional outputs. -------------------------------
1664 RETURN
1665 C-----------------------------------------------------------------------
1666 C Block I.
1667 C The following block handles all error returns due to illegal input
1668 C (ISTATE = -3), as detected before calling the core integrator.
1669 C First the error message routine is called. If the illegal input
1670 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1671 C-----------------------------------------------------------------------
1759 C
1761 RETURN
1762 C
1765 RETURN
1766 C----------------------- End of Subroutine DLSODI ----------------------
1767 END