| blob | 6d8d0a1617f91f49c0b70c036472cf6ad0657002 |
1 *DECK DLSODE
8 C***BEGIN PROLOGUE DLSODE
9 C***PURPOSE Livermore Solver for Ordinary Differential Equations.
10 C DLSODE solves the initial-value problem for stiff or
11 C nonstiff systems of first-order ODE's,
12 C dy/dt = f(t,y), or, in component form,
13 C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)), i=1,...,N.
14 C***CATEGORY I1A
15 C***TYPE DOUBLE PRECISION (SLSODE-S, DLSODE-D)
16 C***KEYWORDS ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM,
17 C STIFF, NONSTIFF
18 C***AUTHOR Hindmarsh, Alan C., (LLNL)
19 C Center for Applied Scientific Computing, L-561
20 C Lawrence Livermore National Laboratory
21 C Livermore, CA 94551.
22 C***DESCRIPTION
23 C
24 C NOTE: The "Usage" and "Arguments" sections treat only a subset of
25 C available options, in condensed fashion. The options
26 C covered and the information supplied will support most
27 C standard uses of DLSODE.
28 C
29 C For more sophisticated uses, full details on all options are
30 C given in the concluding section, headed "Long Description."
31 C A synopsis of the DLSODE Long Description is provided at the
32 C beginning of that section; general topics covered are:
33 C - Elements of the call sequence; optional input and output
34 C - Optional supplemental routines in the DLSODE package
35 C - internal COMMON block
36 C
37 C *Usage:
38 C Communication between the user and the DLSODE package, for normal
39 C situations, is summarized here. This summary describes a subset
40 C of the available options. See "Long Description" for complete
41 C details, including optional communication, nonstandard options,
42 C and instructions for special situations.
43 C
44 C A sample program is given in the "Examples" section.
45 C
46 C Refer to the argument descriptions for the definitions of the
47 C quantities that appear in the following sample declarations.
48 C
49 C For MF = 10,
50 C PARAMETER (LRW = 20 + 16*NEQ, LIW = 20)
51 C For MF = 21 or 22,
52 C PARAMETER (LRW = 22 + 9*NEQ + NEQ**2, LIW = 20 + NEQ)
53 C For MF = 24 or 25,
54 C PARAMETER (LRW = 22 + 10*NEQ + (2*ML+MU)*NEQ,
55 C * LIW = 20 + NEQ)
56 C
57 C EXTERNAL F, JAC
58 C INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW),
59 C * LIW, MF
60 C DOUBLE PRECISION Y(NEQ), T, TOUT, RTOL, ATOL(ntol), RWORK(LRW)
61 C
62 C CALL DLSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
63 C * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
64 C
65 C *Arguments:
66 C F :EXT Name of subroutine for right-hand-side vector f.
67 C This name must be declared EXTERNAL in calling
68 C program. The form of F must be:
69 C
70 C SUBROUTINE F (NEQ, T, Y, YDOT)
71 C INTEGER NEQ
72 C DOUBLE PRECISION T, Y(*), YDOT(*)
73 C
74 C The inputs are NEQ, T, Y. F is to set
75 C
76 C YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)),
77 C i = 1, ..., NEQ .
78 C
79 C NEQ :IN Number of first-order ODE's.
80 C
81 C Y :INOUT Array of values of the y(t) vector, of length NEQ.
82 C Input: For the first call, Y should contain the
83 C values of y(t) at t = T. (Y is an input
84 C variable only if ISTATE = 1.)
85 C Output: On return, Y will contain the values at the
86 C new t-value.
87 C
88 C T :INOUT Value of the independent variable. On return it
89 C will be the current value of t (normally TOUT).
90 C
91 C TOUT :IN Next point where output is desired (.NE. T).
92 C
93 C ITOL :IN 1 or 2 according as ATOL (below) is a scalar or
94 C an array.
95 C
96 C RTOL :IN Relative tolerance parameter (scalar).
97 C
98 C ATOL :IN Absolute tolerance parameter (scalar or array).
99 C If ITOL = 1, ATOL need not be dimensioned.
100 C If ITOL = 2, ATOL must be dimensioned at least NEQ.
101 C
102 C The estimated local error in Y(i) will be controlled
103 C so as to be roughly less (in magnitude) than
104 C
105 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
106 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
107 C
108 C Thus the local error test passes if, in each
109 C component, either the absolute error is less than
110 C ATOL (or ATOL(i)), or the relative error is less
111 C than RTOL.
112 C
113 C Use RTOL = 0.0 for pure absolute error control, and
114 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative
115 C error control. Caution: Actual (global) errors may
116 C exceed these local tolerances, so choose them
117 C conservatively.
118 C
119 C ITASK :IN Flag indicating the task DLSODE is to perform.
120 C Use ITASK = 1 for normal computation of output
121 C values of y at t = TOUT.
122 C
123 C ISTATE:INOUT Index used for input and output to specify the state
124 C of the calculation.
125 C Input:
126 C 1 This is the first call for a problem.
127 C 2 This is a subsequent call.
128 C Output:
129 C 1 Nothing was done, because TOUT was equal to T.
130 C 2 DLSODE was successful (otherwise, negative).
131 C Note that ISTATE need not be modified after a
132 C successful return.
133 C -1 Excess work done on this call (perhaps wrong
134 C MF).
135 C -2 Excess accuracy requested (tolerances too
136 C small).
137 C -3 Illegal input detected (see printed message).
138 C -4 Repeated error test failures (check all
139 C inputs).
140 C -5 Repeated convergence failures (perhaps bad
141 C Jacobian supplied or wrong choice of MF or
142 C tolerances).
143 C -6 Error weight became zero during problem
144 C (solution component i vanished, and ATOL or
145 C ATOL(i) = 0.).
146 C
147 C IOPT :IN Flag indicating whether optional inputs are used:
148 C 0 No.
149 C 1 Yes. (See "Optional inputs" under "Long
150 C Description," Part 1.)
151 C
152 C RWORK :WORK Real work array of length at least:
153 C 20 + 16*NEQ for MF = 10,
154 C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
155 C 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25.
156 C
157 C LRW :IN Declared length of RWORK (in user's DIMENSION
158 C statement).
159 C
160 C IWORK :WORK Integer work array of length at least:
161 C 20 for MF = 10,
162 C 20 + NEQ for MF = 21, 22, 24, or 25.
163 C
164 C If MF = 24 or 25, input in IWORK(1),IWORK(2) the
165 C lower and upper Jacobian half-bandwidths ML,MU.
166 C
167 C On return, IWORK contains information that may be
168 C of interest to the user:
169 C
170 C Name Location Meaning
171 C ----- --------- -----------------------------------------
172 C NST IWORK(11) Number of steps taken for the problem so
173 C far.
174 C NFE IWORK(12) Number of f evaluations for the problem
175 C so far.
176 C NJE IWORK(13) Number of Jacobian evaluations (and of
177 C matrix LU decompositions) for the problem
178 C so far.
179 C NQU IWORK(14) Method order last used (successfully).
180 C LENRW IWORK(17) Length of RWORK actually required. This
181 C is defined on normal returns and on an
182 C illegal input return for insufficient
183 C storage.
184 C LENIW IWORK(18) Length of IWORK actually required. This
185 C is defined on normal returns and on an
186 C illegal input return for insufficient
187 C storage.
188 C
189 C LIW :IN Declared length of IWORK (in user's DIMENSION
190 C statement).
191 C
192 C JAC :EXT Name of subroutine for Jacobian matrix (MF =
193 C 21 or 24). If used, this name must be declared
194 C EXTERNAL in calling program. If not used, pass a
195 C dummy name. The form of JAC must be:
196 C
197 C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
198 C INTEGER NEQ, ML, MU, NROWPD
199 C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
200 C
201 C See item c, under "Description" below for more
202 C information about JAC.
203 C
204 C MF :IN Method flag. Standard values are:
205 C 10 Nonstiff (Adams) method, no Jacobian used.
206 C 21 Stiff (BDF) method, user-supplied full Jacobian.
207 C 22 Stiff method, internally generated full
208 C Jacobian.
209 C 24 Stiff method, user-supplied banded Jacobian.
210 C 25 Stiff method, internally generated banded
211 C Jacobian.
212 C
213 C *Description:
214 C DLSODE solves the initial value problem for stiff or nonstiff
215 C systems of first-order ODE's,
216 C
217 C dy/dt = f(t,y) ,
218 C
219 C or, in component form,
220 C
221 C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ))
222 C (i = 1, ..., NEQ) .
223 C
224 C DLSODE is a package based on the GEAR and GEARB packages, and on
225 C the October 23, 1978, version of the tentative ODEPACK user
226 C interface standard, with minor modifications.
227 C
228 C The steps in solving such a problem are as follows.
229 C
230 C a. First write a subroutine of the form
231 C
232 C SUBROUTINE F (NEQ, T, Y, YDOT)
233 C INTEGER NEQ
234 C DOUBLE PRECISION T, Y(*), YDOT(*)
235 C
236 C which supplies the vector function f by loading YDOT(i) with
237 C f(i).
238 C
239 C b. Next determine (or guess) whether or not the problem is stiff.
240 C Stiffness occurs when the Jacobian matrix df/dy has an
241 C eigenvalue whose real part is negative and large in magnitude
242 C compared to the reciprocal of the t span of interest. If the
243 C problem is nonstiff, use method flag MF = 10. If it is stiff,
244 C there are four standard choices for MF, and DLSODE requires the
245 C Jacobian matrix in some form. This matrix is regarded either
246 C as full (MF = 21 or 22), or banded (MF = 24 or 25). In the
247 C banded case, DLSODE requires two half-bandwidth parameters ML
248 C and MU. These are, respectively, the widths of the lower and
249 C upper parts of the band, excluding the main diagonal. Thus the
250 C band consists of the locations (i,j) with
251 C
252 C i - ML <= j <= i + MU ,
253 C
254 C and the full bandwidth is ML + MU + 1 .
255 C
256 C c. If the problem is stiff, you are encouraged to supply the
257 C Jacobian directly (MF = 21 or 24), but if this is not feasible,
258 C DLSODE will compute it internally by difference quotients (MF =
259 C 22 or 25). If you are supplying the Jacobian, write a
260 C subroutine of the form
261 C
262 C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
263 C INTEGER NEQ, ML, MU, NRWOPD
264 C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
265 C
266 C which provides df/dy by loading PD as follows:
267 C - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j),
268 C the partial derivative of f(i) with respect to y(j). (Ignore
269 C the ML and MU arguments in this case.)
270 C - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with
271 C df(i)/dy(j); i.e., load the diagonal lines of df/dy into the
272 C rows of PD from the top down.
273 C - In either case, only nonzero elements need be loaded.
274 C
275 C d. Write a main program that calls subroutine DLSODE once for each
276 C point at which answers are desired. This should also provide
277 C for possible use of logical unit 6 for output of error messages
278 C by DLSODE.
279 C
280 C Before the first call to DLSODE, set ISTATE = 1, set Y and T to
281 C the initial values, and set TOUT to the first output point. To
282 C continue the integration after a successful return, simply
283 C reset TOUT and call DLSODE again. No other parameters need be
284 C reset.
285 C
286 C *Examples:
287 C The following is a simple example problem, with the coding needed
288 C for its solution by DLSODE. The problem is from chemical kinetics,
289 C and consists of the following three rate equations:
290 C
291 C dy1/dt = -.04*y1 + 1.E4*y2*y3
292 C dy2/dt = .04*y1 - 1.E4*y2*y3 - 3.E7*y2**2
293 C dy3/dt = 3.E7*y2**2
294 C
295 C on the interval from t = 0.0 to t = 4.E10, with initial conditions
296 C y1 = 1.0, y2 = y3 = 0. The problem is stiff.
297 C
298 C The following coding solves this problem with DLSODE, using
299 C MF = 21 and printing results at t = .4, 4., ..., 4.E10. It uses
300 C ITOL = 2 and ATOL much smaller for y2 than for y1 or y3 because y2
301 C has much smaller values. At the end of the run, statistical
302 C quantities of interest are printed.
303 C
304 C EXTERNAL FEX, JEX
305 C INTEGER IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW,
306 C * MF, NEQ
307 C DOUBLE PRECISION ATOL(3), RTOL, RWORK(58), T, TOUT, Y(3)
308 C NEQ = 3
309 C Y(1) = 1.D0
310 C Y(2) = 0.D0
311 C Y(3) = 0.D0
312 C T = 0.D0
313 C TOUT = .4D0
314 C ITOL = 2
315 C RTOL = 1.D-4
316 C ATOL(1) = 1.D-6
317 C ATOL(2) = 1.D-10
318 C ATOL(3) = 1.D-6
319 C ITASK = 1
320 C ISTATE = 1
321 C IOPT = 0
322 C LRW = 58
323 C LIW = 23
324 C MF = 21
325 C DO 40 IOUT = 1,12
326 C CALL DLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
327 C * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
328 C WRITE(6,20) T, Y(1), Y(2), Y(3)
329 C 20 FORMAT(' At t =',D12.4,' y =',3D14.6)
330 C IF (ISTATE .LT. 0) GO TO 80
331 C 40 TOUT = TOUT*10.D0
332 C WRITE(6,60) IWORK(11), IWORK(12), IWORK(13)
333 C 60 FORMAT(/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)
334 C STOP
335 C 80 WRITE(6,90) ISTATE
336 C 90 FORMAT(///' Error halt.. ISTATE =',I3)
337 C STOP
338 C END
339 C
340 C SUBROUTINE FEX (NEQ, T, Y, YDOT)
341 C INTEGER NEQ
342 C DOUBLE PRECISION T, Y(3), YDOT(3)
343 C YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
344 C YDOT(3) = 3.D7*Y(2)*Y(2)
345 C YDOT(2) = -YDOT(1) - YDOT(3)
346 C RETURN
347 C END
348 C
349 C SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD)
350 C INTEGER NEQ, ML, MU, NRPD
351 C DOUBLE PRECISION T, Y(3), PD(NRPD,3)
352 C PD(1,1) = -.04D0
353 C PD(1,2) = 1.D4*Y(3)
354 C PD(1,3) = 1.D4*Y(2)
355 C PD(2,1) = .04D0
356 C PD(2,3) = -PD(1,3)
357 C PD(3,2) = 6.D7*Y(2)
358 C PD(2,2) = -PD(1,2) - PD(3,2)
359 C RETURN
360 C END
361 C
362 C The output from this program (on a Cray-1 in single precision)
363 C is as follows.
364 C
365 C At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02
366 C At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02
367 C At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01
368 C At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01
369 C At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01
370 C At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01
371 C At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01
372 C At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01
373 C At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01
374 C At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01
375 C At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01
376 C At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00
377 C
378 C No. steps = 330, No. f-s = 405, No. J-s = 69
379 C
380 C *Accuracy:
381 C The accuracy of the solution depends on the choice of tolerances
382 C RTOL and ATOL. Actual (global) errors may exceed these local
383 C tolerances, so choose them conservatively.
384 C
385 C *Cautions:
386 C The work arrays should not be altered between calls to DLSODE for
387 C the same problem, except possibly for the conditional and optional
388 C inputs.
389 C
390 C *Portability:
391 C Since NEQ is dimensioned inside DLSODE, some compilers may object
392 C to a call to DLSODE with NEQ a scalar variable. In this event,
393 C use DIMENSION NEQ(1). Similar remarks apply to RTOL and ATOL.
394 C
395 C Note to Cray users:
396 C For maximum efficiency, use the CFT77 compiler. Appropriate
397 C compiler optimization directives have been inserted for CFT77.
398 C
399 C *Reference:
400 C Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
401 C Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds.
402 C (North-Holland, Amsterdam, 1983), pp. 55-64.
403 C
404 C *Long Description:
405 C The following complete description of the user interface to
406 C DLSODE consists of four parts:
407 C
408 C 1. The call sequence to subroutine DLSODE, which is a driver
409 C routine for the solver. This includes descriptions of both
410 C the call sequence arguments and user-supplied routines.
411 C Following these descriptions is a description of optional
412 C inputs available through the call sequence, and then a
413 C description of optional outputs in the work arrays.
414 C
415 C 2. Descriptions of other routines in the DLSODE package that may
416 C be (optionally) called by the user. These provide the ability
417 C to alter error message handling, save and restore the internal
418 C COMMON, and obtain specified derivatives of the solution y(t).
419 C
420 C 3. Descriptions of COMMON block to be declared in overlay or
421 C similar environments, or to be saved when doing an interrupt
422 C of the problem and continued solution later.
423 C
424 C 4. Description of two routines in the DLSODE package, either of
425 C which the user may replace with his own version, if desired.
426 C These relate to the measurement of errors.
427 C
428 C
429 C Part 1. Call Sequence
430 C ----------------------
431 C
432 C Arguments
433 C ---------
434 C The call sequence parameters used for input only are
435 C
436 C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
437 C
438 C and those used for both input and output are
439 C
440 C Y, T, ISTATE.
441 C
442 C The work arrays RWORK and IWORK are also used for conditional and
443 C optional inputs and optional outputs. (The term output here
444 C refers to the return from subroutine DLSODE to the user's calling
445 C program.)
446 C
447 C The legality of input parameters will be thoroughly checked on the
448 C initial call for the problem, but not checked thereafter unless a
449 C change in input parameters is flagged by ISTATE = 3 on input.
450 C
451 C The descriptions of the call arguments are as follows.
452 C
453 C F The name of the user-supplied subroutine defining the ODE
454 C system. The system must be put in the first-order form
455 C dy/dt = f(t,y), where f is a vector-valued function of
456 C the scalar t and the vector y. Subroutine F is to compute
457 C the function f. It is to have the form
458 C
459 C SUBROUTINE F (NEQ, T, Y, YDOT)
460 C DOUBLE PRECISION T, Y(*), YDOT(*)
461 C
462 C where NEQ, T, and Y are input, and the array YDOT =
463 C f(T,Y) is output. Y and YDOT are arrays of length NEQ.
464 C Subroutine F should not alter Y(1),...,Y(NEQ). F must be
465 C declared EXTERNAL in the calling program.
466 C
467 C Subroutine F may access user-defined quantities in
468 C NEQ(2),... and/or in Y(NEQ(1)+1),..., if NEQ is an array
469 C (dimensioned in F) and/or Y has length exceeding NEQ(1).
470 C See the descriptions of NEQ and Y below.
471 C
472 C If quantities computed in the F routine are needed
473 C externally to DLSODE, an extra call to F should be made
474 C for this purpose, for consistent and accurate results.
475 C If only the derivative dy/dt is needed, use DINTDY
476 C instead.
477 C
478 C NEQ The size of the ODE system (number of first-order
479 C ordinary differential equations). Used only for input.
480 C NEQ may be decreased, but not increased, during the
481 C problem. If NEQ is decreased (with ISTATE = 3 on input),
482 C the remaining components of Y should be left undisturbed,
483 C if these are to be accessed in F and/or JAC.
484 C
485 C Normally, NEQ is a scalar, and it is generally referred
486 C to as a scalar in this user interface description.
487 C However, NEQ may be an array, with NEQ(1) set to the
488 C system size. (The DLSODE package accesses only NEQ(1).)
489 C In either case, this parameter is passed as the NEQ
490 C argument in all calls to F and JAC. Hence, if it is an
491 C array, locations NEQ(2),... may be used to store other
492 C integer data and pass it to F and/or JAC. Subroutines
493 C F and/or JAC must include NEQ in a DIMENSION statement
494 C in that case.
495 C
496 C Y A real array for the vector of dependent variables, of
497 C length NEQ or more. Used for both input and output on
498 C the first call (ISTATE = 1), and only for output on
499 C other calls. On the first call, Y must contain the
500 C vector of initial values. On output, Y contains the
501 C computed solution vector, evaluated at T. If desired,
502 C the Y array may be used for other purposes between
503 C calls to the solver.
504 C
505 C This array is passed as the Y argument in all calls to F
506 C and JAC. Hence its length may exceed NEQ, and locations
507 C Y(NEQ+1),... may be used to store other real data and
508 C pass it to F and/or JAC. (The DLSODE package accesses
509 C only Y(1),...,Y(NEQ).)
510 C
511 C T The independent variable. On input, T is used only on
512 C the first call, as the initial point of the integration.
513 C On output, after each call, T is the value at which a
514 C computed solution Y is evaluated (usually the same as
515 C TOUT). On an error return, T is the farthest point
516 C reached.
517 C
518 C TOUT The next value of T at which a computed solution is
519 C desired. Used only for input.
520 C
521 C When starting the problem (ISTATE = 1), TOUT may be equal
522 C to T for one call, then should not equal T for the next
523 C call. For the initial T, an input value of TOUT .NE. T
524 C is used in order to determine the direction of the
525 C integration (i.e., the algebraic sign of the step sizes)
526 C and the rough scale of the problem. Integration in
527 C either direction (forward or backward in T) is permitted.
528 C
529 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored
530 C after the first call (i.e., the first call with
531 C TOUT .NE. T). Otherwise, TOUT is required on every call.
532 C
533 C If ITASK = 1, 3, or 4, the values of TOUT need not be
534 C monotone, but a value of TOUT which backs up is limited
535 C to the current internal T interval, whose endpoints are
536 C TCUR - HU and TCUR. (See "Optional Outputs" below for
537 C TCUR and HU.)
538 C
539 C
540 C ITOL An indicator for the type of error control. See
541 C description below under ATOL. Used only for input.
542 C
543 C RTOL A relative error tolerance parameter, either a scalar or
544 C an array of length NEQ. See description below under
545 C ATOL. Input only.
546 C
547 C ATOL An absolute error tolerance parameter, either a scalar or
548 C an array of length NEQ. Input only.
549 C
550 C The input parameters ITOL, RTOL, and ATOL determine the
551 C error control performed by the solver. The solver will
552 C control the vector e = (e(i)) of estimated local errors
553 C in Y, according to an inequality of the form
554 C
555 C rms-norm of ( e(i)/EWT(i) ) <= 1,
556 C
557 C where
558 C
559 C EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
560 C
561 C and the rms-norm (root-mean-square norm) here is
562 C
563 C rms-norm(v) = SQRT(sum v(i)**2 / NEQ).
564 C
565 C Here EWT = (EWT(i)) is a vector of weights which must
566 C always be positive, and the values of RTOL and ATOL
567 C should all be nonnegative. The following table gives the
568 C types (scalar/array) of RTOL and ATOL, and the
569 C corresponding form of EWT(i).
570 C
571 C ITOL RTOL ATOL EWT(i)
572 C ---- ------ ------ -----------------------------
573 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
574 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
575 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
576 C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i)
577 C
578 C When either of these parameters is a scalar, it need not
579 C be dimensioned in the user's calling program.
580 C
581 C If none of the above choices (with ITOL, RTOL, and ATOL
582 C fixed throughout the problem) is suitable, more general
583 C error controls can be obtained by substituting
584 C user-supplied routines for the setting of EWT and/or for
585 C the norm calculation. See Part 4 below.
586 C
587 C If global errors are to be estimated by making a repeated
588 C run on the same problem with smaller tolerances, then all
589 C components of RTOL and ATOL (i.e., of EWT) should be
590 C scaled down uniformly.
591 C
592 C ITASK An index specifying the task to be performed. Input
593 C only. ITASK has the following values and meanings:
594 C 1 Normal computation of output values of y(t) at
595 C t = TOUT (by overshooting and interpolating).
596 C 2 Take one step only and return.
597 C 3 Stop at the first internal mesh point at or beyond
598 C t = TOUT and return.
599 C 4 Normal computation of output values of y(t) at
600 C t = TOUT but without overshooting t = TCRIT. TCRIT
601 C must be input as RWORK(1). TCRIT may be equal to or
602 C beyond TOUT, but not behind it in the direction of
603 C integration. This option is useful if the problem
604 C has a singularity at or beyond t = TCRIT.
605 C 5 Take one step, without passing TCRIT, and return.
606 C TCRIT must be input as RWORK(1).
607 C
608 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
609 C (within roundoff), it will return T = TCRIT (exactly) to
610 C indicate this (unless ITASK = 4 and TOUT comes before
611 C TCRIT, in which case answers at T = TOUT are returned
612 C first).
613 C
614 C ISTATE An index used for input and output to specify the state
615 C of the calculation.
616 C
617 C On input, the values of ISTATE are as follows:
618 C 1 This is the first call for the problem
619 C (initializations will be done). See "Note" below.
620 C 2 This is not the first call, and the calculation is to
621 C continue normally, with no change in any input
622 C parameters except possibly TOUT and ITASK. (If ITOL,
623 C RTOL, and/or ATOL are changed between calls with
624 C ISTATE = 2, the new values will be used but not
625 C tested for legality.)
626 C 3 This is not the first call, and the calculation is to
627 C continue normally, but with a change in input
628 C parameters other than TOUT and ITASK. Changes are
629 C allowed in NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF,
630 C ML, MU, and any of the optional inputs except H0.
631 C (See IWORK description for ML and MU.)
632 C
633 C Note: A preliminary call with TOUT = T is not counted as
634 C a first call here, as no initialization or checking of
635 C input is done. (Such a call is sometimes useful for the
636 C purpose of outputting the initial conditions.) Thus the
637 C first call for which TOUT .NE. T requires ISTATE = 1 on
638 C input.
639 C
640 C On output, ISTATE has the following values and meanings:
641 C 1 Nothing was done, as TOUT was equal to T with
642 C ISTATE = 1 on input.
643 C 2 The integration was performed successfully.
644 C -1 An excessive amount of work (more than MXSTEP steps)
645 C was done on this call, before completing the
646 C requested task, but the integration was otherwise
647 C successful as far as T. (MXSTEP is an optional input
648 C and is normally 500.) To continue, the user may
649 C simply reset ISTATE to a value >1 and call again (the
650 C excess work step counter will be reset to 0). In
651 C addition, the user may increase MXSTEP to avoid this
652 C error return; see "Optional Inputs" below.
653 C -2 Too much accuracy was requested for the precision of
654 C the machine being used. This was detected before
655 C completing the requested task, but the integration
656 C was successful as far as T. To continue, the
657 C tolerance parameters must be reset, and ISTATE must
658 C be set to 3. The optional output TOLSF may be used
659 C for this purpose. (Note: If this condition is
660 C detected before taking any steps, then an illegal
661 C input return (ISTATE = -3) occurs instead.)
662 C -3 Illegal input was detected, before taking any
663 C integration steps. See written message for details.
664 C (Note: If the solver detects an infinite loop of
665 C calls to the solver with illegal input, it will cause
666 C the run to stop.)
667 C -4 There were repeated error-test failures on one
668 C attempted step, before completing the requested task,
669 C but the integration was successful as far as T. The
670 C problem may have a singularity, or the input may be
671 C inappropriate.
672 C -5 There were repeated convergence-test failures on one
673 C attempted step, before completing the requested task,
674 C but the integration was successful as far as T. This
675 C may be caused by an inaccurate Jacobian matrix, if
676 C one is being used.
677 C -6 EWT(i) became zero for some i during the integration.
678 C Pure relative error control (ATOL(i)=0.0) was
679 C requested on a variable which has now vanished. The
680 C integration was successful as far as T.
681 C
682 C Note: Since the normal output value of ISTATE is 2, it
683 C does not need to be reset for normal continuation. Also,
684 C since a negative input value of ISTATE will be regarded
685 C as illegal, a negative output value requires the user to
686 C change it, and possibly other inputs, before calling the
687 C solver again.
688 C
689 C IOPT An integer flag to specify whether any optional inputs
690 C are being used on this call. Input only. The optional
691 C inputs are listed under a separate heading below.
692 C 0 No optional inputs are being used. Default values
693 C will be used in all cases.
694 C 1 One or more optional inputs are being used.
695 C
696 C RWORK A real working array (double precision). The length of
697 C RWORK must be at least
698 C
699 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM
700 C
701 C where
702 C NYH = the initial value of NEQ,
703 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
704 C smaller value is given as an optional input),
705 C LWM = 0 if MITER = 0,
706 C LWM = NEQ**2 + 2 if MITER = 1 or 2,
707 C LWM = NEQ + 2 if MITER = 3, and
708 C LWM = (2*ML + MU + 1)*NEQ + 2
709 C if MITER = 4 or 5.
710 C (See the MF description below for METH and MITER.)
711 C
712 C Thus if MAXORD has its default value and NEQ is constant,
713 C this length is:
714 C 20 + 16*NEQ for MF = 10,
715 C 22 + 16*NEQ + NEQ**2 for MF = 11 or 12,
716 C 22 + 17*NEQ for MF = 13,
717 C 22 + 17*NEQ + (2*ML + MU)*NEQ for MF = 14 or 15,
718 C 20 + 9*NEQ for MF = 20,
719 C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
720 C 22 + 10*NEQ for MF = 23,
721 C 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25.
722 C
723 C The first 20 words of RWORK are reserved for conditional
724 C and optional inputs and optional outputs.
725 C
726 C The following word in RWORK is a conditional input:
727 C RWORK(1) = TCRIT, the critical value of t which the
728 C solver is not to overshoot. Required if ITASK
729 C is 4 or 5, and ignored otherwise. See ITASK.
730 C
731 C LRW The length of the array RWORK, as declared by the user.
732 C (This will be checked by the solver.)
733 C
734 C IWORK An integer work array. Its length must be at least
735 C 20 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or
736 C 20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25).
737 C (See the MF description below for MITER.) The first few
738 C words of IWORK are used for conditional and optional
739 C inputs and optional outputs.
740 C
741 C The following two words in IWORK are conditional inputs:
742 C IWORK(1) = ML These are the lower and upper half-
743 C IWORK(2) = MU bandwidths, respectively, of the banded
744 C Jacobian, excluding the main diagonal.
745 C The band is defined by the matrix locations
746 C (i,j) with i - ML <= j <= i + MU. ML and MU
747 C must satisfy 0 <= ML,MU <= NEQ - 1. These are
748 C required if MITER is 4 or 5, and ignored
749 C otherwise. ML and MU may in fact be the band
750 C parameters for a matrix to which df/dy is only
751 C approximately equal.
752 C
753 C LIW The length of the array IWORK, as declared by the user.
754 C (This will be checked by the solver.)
755 C
756 C Note: The work arrays must not be altered between calls to DLSODE
757 C for the same problem, except possibly for the conditional and
758 C optional inputs, and except for the last 3*NEQ words of RWORK.
759 C The latter space is used for internal scratch space, and so is
760 C available for use by the user outside DLSODE between calls, if
761 C desired (but not for use by F or JAC).
762 C
763 C JAC The name of the user-supplied routine (MITER = 1 or 4) to
764 C compute the Jacobian matrix, df/dy, as a function of the
765 C scalar t and the vector y. (See the MF description below
766 C for MITER.) It is to have the form
767 C
768 C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
769 C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
770 C
771 C where NEQ, T, Y, ML, MU, and NROWPD are input and the
772 C array PD is to be loaded with partial derivatives
773 C (elements of the Jacobian matrix) on output. PD must be
774 C given a first dimension of NROWPD. T and Y have the same
775 C meaning as in subroutine F.
776 C
777 C In the full matrix case (MITER = 1), ML and MU are
778 C ignored, and the Jacobian is to be loaded into PD in
779 C columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
780 C
781 C In the band matrix case (MITER = 4), the elements within
782 C the band are to be loaded into PD in columnwise manner,
783 C with diagonal lines of df/dy loaded into the rows of PD.
784 C Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML
785 C and MU are the half-bandwidth parameters (see IWORK).
786 C The locations in PD in the two triangular areas which
787 C correspond to nonexistent matrix elements can be ignored
788 C or loaded arbitrarily, as they are overwritten by DLSODE.
789 C
790 C JAC need not provide df/dy exactly. A crude approximation
791 C (possibly with a smaller bandwidth) will do.
792 C
793 C In either case, PD is preset to zero by the solver, so
794 C that only the nonzero elements need be loaded by JAC.
795 C Each call to JAC is preceded by a call to F with the same
796 C arguments NEQ, T, and Y. Thus to gain some efficiency,
797 C intermediate quantities shared by both calculations may
798 C be saved in a user COMMON block by F and not recomputed
799 C by JAC, if desired. Also, JAC may alter the Y array, if
800 C desired. JAC must be declared EXTERNAL in the calling
801 C program.
802 C
803 C Subroutine JAC may access user-defined quantities in
804 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
805 C (dimensioned in JAC) and/or Y has length exceeding
806 C NEQ(1). See the descriptions of NEQ and Y above.
807 C
808 C MF The method flag. Used only for input. The legal values
809 C of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24,
810 C and 25. MF has decimal digits METH and MITER:
811 C MF = 10*METH + MITER .
812 C
813 C METH indicates the basic linear multistep method:
814 C 1 Implicit Adams method.
815 C 2 Method based on backward differentiation formulas
816 C (BDF's).
817 C
818 C MITER indicates the corrector iteration method:
819 C 0 Functional iteration (no Jacobian matrix is
820 C involved).
821 C 1 Chord iteration with a user-supplied full (NEQ by
822 C NEQ) Jacobian.
823 C 2 Chord iteration with an internally generated
824 C (difference quotient) full Jacobian (using NEQ
825 C extra calls to F per df/dy value).
826 C 3 Chord iteration with an internally generated
827 C diagonal Jacobian approximation (using one extra call
828 C to F per df/dy evaluation).
829 C 4 Chord iteration with a user-supplied banded Jacobian.
830 C 5 Chord iteration with an internally generated banded
831 C Jacobian (using ML + MU + 1 extra calls to F per
832 C df/dy evaluation).
833 C
834 C If MITER = 1 or 4, the user must supply a subroutine JAC
835 C (the name is arbitrary) as described above under JAC.
836 C For other values of MITER, a dummy argument can be used.
837 C
838 C Optional Inputs
839 C ---------------
840 C The following is a list of the optional inputs provided for in the
841 C call sequence. (See also Part 2.) For each such input variable,
842 C this table lists its name as used in this documentation, its
843 C location in the call sequence, its meaning, and the default value.
844 C The use of any of these inputs requires IOPT = 1, and in that case
845 C all of these inputs are examined. A value of zero for any of
846 C these optional inputs will cause the default value to be used.
847 C Thus to use a subset of the optional inputs, simply preload
848 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively,
849 C and then set those of interest to nonzero values.
850 C
851 C Name Location Meaning and default value
852 C ------ --------- -----------------------------------------------
853 C H0 RWORK(5) Step size to be attempted on the first step.
854 C The default value is determined by the solver.
855 C HMAX RWORK(6) Maximum absolute step size allowed. The
856 C default value is infinite.
857 C HMIN RWORK(7) Minimum absolute step size allowed. The
858 C default value is 0. (This lower bound is not
859 C enforced on the final step before reaching
860 C TCRIT when ITASK = 4 or 5.)
861 C MAXORD IWORK(5) Maximum order to be allowed. The default value
862 C is 12 if METH = 1, and 5 if METH = 2. (See the
863 C MF description above for METH.) If MAXORD
864 C exceeds the default value, it will be reduced
865 C to the default value. If MAXORD is changed
866 C during the problem, it may cause the current
867 C order to be reduced.
868 C MXSTEP IWORK(6) Maximum number of (internally defined) steps
869 C allowed during one call to the solver. The
870 C default value is 500.
871 C MXHNIL IWORK(7) Maximum number of messages printed (per
872 C problem) warning that T + H = T on a step
873 C (H = step size). This must be positive to
874 C result in a nondefault value. The default
875 C value is 10.
876 C
877 C Optional Outputs
878 C ----------------
879 C As optional additional output from DLSODE, the variables listed
880 C below are quantities related to the performance of DLSODE which
881 C are available to the user. These are communicated by way of the
882 C work arrays, but also have internal mnemonic names as shown.
883 C Except where stated otherwise, all of these outputs are defined on
884 C any successful return from DLSODE, and on any return with ISTATE =
885 C -1, -2, -4, -5, or -6. On an illegal input return (ISTATE = -3),
886 C they will be unchanged from their existing values (if any), except
887 C possibly for TOLSF, LENRW, and LENIW. On any error return,
888 C outputs relevant to the error will be defined, as noted below.
889 C
890 C Name Location Meaning
891 C ----- --------- ------------------------------------------------
892 C HU RWORK(11) Step size in t last used (successfully).
893 C HCUR RWORK(12) Step size to be attempted on the next step.
894 C TCUR RWORK(13) Current value of the independent variable which
895 C the solver has actually reached, i.e., the
896 C current internal mesh point in t. On output,
897 C TCUR will always be at least as far as the
898 C argument T, but may be farther (if interpolation
899 C was done).
900 C TOLSF RWORK(14) Tolerance scale factor, greater than 1.0,
901 C computed when a request for too much accuracy
902 C was detected (ISTATE = -3 if detected at the
903 C start of the problem, ISTATE = -2 otherwise).
904 C If ITOL is left unaltered but RTOL and ATOL are
905 C uniformly scaled up by a factor of TOLSF for the
906 C next call, then the solver is deemed likely to
907 C succeed. (The user may also ignore TOLSF and
908 C alter the tolerance parameters in any other way
909 C appropriate.)
910 C NST IWORK(11) Number of steps taken for the problem so far.
911 C NFE IWORK(12) Number of F evaluations for the problem so far.
912 C NJE IWORK(13) Number of Jacobian evaluations (and of matrix LU
913 C decompositions) for the problem so far.
914 C NQU IWORK(14) Method order last used (successfully).
915 C NQCUR IWORK(15) Order to be attempted on the next step.
916 C IMXER IWORK(16) Index of the component of largest magnitude in
917 C the weighted local error vector ( e(i)/EWT(i) ),
918 C on an error return with ISTATE = -4 or -5.
919 C LENRW IWORK(17) Length of RWORK actually required. This is
920 C defined on normal returns and on an illegal
921 C input return for insufficient storage.
922 C LENIW IWORK(18) Length of IWORK actually required. This is
923 C defined on normal returns and on an illegal
924 C input return for insufficient storage.
925 C
926 C The following two arrays are segments of the RWORK array which may
927 C also be of interest to the user as optional outputs. For each
928 C array, the table below gives its internal name, its base address
929 C in RWORK, and its description.
930 C
931 C Name Base address Description
932 C ---- ------------ ----------------------------------------------
933 C YH 21 The Nordsieck history array, of size NYH by
934 C (NQCUR + 1), where NYH is the initial value of
935 C NEQ. For j = 0,1,...,NQCUR, column j + 1 of
936 C YH contains HCUR**j/factorial(j) times the jth
937 C derivative of the interpolating polynomial
938 C currently representing the solution, evaluated
939 C at t = TCUR.
940 C ACOR LENRW-NEQ+1 Array of size NEQ used for the accumulated
941 C corrections on each step, scaled on output to
942 C represent the estimated local error in Y on
943 C the last step. This is the vector e in the
944 C description of the error control. It is
945 C defined only on successful return from DLSODE.
946 C
947 C
948 C Part 2. Other Callable Routines
949 C --------------------------------
950 C
951 C The following are optional calls which the user may make to gain
952 C additional capabilities in conjunction with DLSODE.
953 C
954 C Form of call Function
955 C ------------------------ ----------------------------------------
956 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
957 C output of messages from DLSODE, if the
958 C default is not desired. The default
959 C value of LUN is 6. This call may be made
960 C at any time and will take effect
961 C immediately.
962 C CALL XSETF(MFLAG) Set a flag to control the printing of
963 C messages by DLSODE. MFLAG = 0 means do
964 C not print. (Danger: this risks losing
965 C valuable information.) MFLAG = 1 means
966 C print (the default). This call may be
967 C made at any time and will take effect
968 C immediately.
969 C CALL DSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the
970 C internal COMMON blocks used by DLSODE
971 C (see Part 3 below). RSAV must be a
972 C real array of length 218 or more, and
973 C ISAV must be an integer array of length
974 C 37 or more. JOB = 1 means save COMMON
975 C into RSAV/ISAV. JOB = 2 means restore
976 C COMMON from same. DSRCOM is useful if
977 C one is interrupting a run and restarting
978 C later, or alternating between two or
979 C more problems solved with DLSODE.
980 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
981 C (see below) orders, at a specified point t, if
982 C desired. It may be called only after a
983 C successful return from DLSODE. Detailed
984 C instructions follow.
985 C
986 C Detailed instructions for using DINTDY
987 C --------------------------------------
988 C The form of the CALL is:
989 C
990 C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
991 C
992 C The input parameters are:
993 C
994 C T Value of independent variable where answers are
995 C desired (normally the same as the T last returned by
996 C DLSODE). For valid results, T must lie between
997 C TCUR - HU and TCUR. (See "Optional Outputs" above
998 C for TCUR and HU.)
999 C K Integer order of the derivative desired. K must
1000 C satisfy 0 <= K <= NQCUR, where NQCUR is the current
1001 C order (see "Optional Outputs"). The capability
1002 C corresponding to K = 0, i.e., computing y(t), is
1003 C already provided by DLSODE directly. Since
1004 C NQCUR >= 1, the first derivative dy/dt is always
1005 C available with DINTDY.
1006 C RWORK(21) The base address of the history array YH.
1007 C NYH Column length of YH, equal to the initial value of NEQ.
1008 C
1009 C The output parameters are:
1010 C
1011 C DKY Real array of length NEQ containing the computed value
1012 C of the Kth derivative of y(t).
1013 C IFLAG Integer flag, returned as 0 if K and T were legal,
1014 C -1 if K was illegal, and -2 if T was illegal.
1015 C On an error return, a message is also written.
1016 C
1017 C
1018 C Part 3. Common Blocks
1019 C ----------------------
1020 C
1021 C If DLSODE is to be used in an overlay situation, the user must
1022 C declare, in the primary overlay, the variables in:
1023 C (1) the call sequence to DLSODE,
1024 C (2) the internal COMMON block /DLS001/, of length 255
1025 C (218 double precision words followed by 37 integer words).
1026 C
1027 C If DLSODE is used on a system in which the contents of internal
1028 C COMMON blocks are not preserved between calls, the user should
1029 C declare the above COMMON block in his main program to insure that
1030 C its contents are preserved.
1031 C
1032 C If the solution of a given problem by DLSODE is to be interrupted
1033 C and then later continued, as when restarting an interrupted run or
1034 C alternating between two or more problems, the user should save,
1035 C following the return from the last DLSODE call prior to the
1036 C interruption, the contents of the call sequence variables and the
1037 C internal COMMON block, and later restore these values before the
1038 C next DLSODE call for that problem. In addition, if XSETUN and/or
1039 C XSETF was called for non-default handling of error messages, then
1040 C these calls must be repeated. To save and restore the COMMON
1041 C block, use subroutine DSRCOM (see Part 2 above).
1042 C
1043 C
1044 C Part 4. Optionally Replaceable Solver Routines
1045 C -----------------------------------------------
1046 C
1047 C Below are descriptions of two routines in the DLSODE package which
1048 C relate to the measurement of errors. Either routine can be
1049 C replaced by a user-supplied version, if desired. However, since
1050 C such a replacement may have a major impact on performance, it
1051 C should be done only when absolutely necessary, and only with great
1052 C caution. (Note: The means by which the package version of a
1053 C routine is superseded by the user's version may be system-
1054 C dependent.)
1055 C
1056 C DEWSET
1057 C ------
1058 C The following subroutine is called just before each internal
1059 C integration step, and sets the array of error weights, EWT, as
1060 C described under ITOL/RTOL/ATOL above:
1061 C
1062 C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1063 C
1064 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODE call
1065 C sequence, YCUR contains the current dependent variable vector,
1066 C and EWT is the array of weights set by DEWSET.
1067 C
1068 C If the user supplies this subroutine, it must return in EWT(i)
1069 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1070 C in Y(i) to. The EWT array returned by DEWSET is passed to the
1071 C DVNORM routine (see below), and also used by DLSODE in the
1072 C computation of the optional output IMXER, the diagonal Jacobian
1073 C approximation, and the increments for difference quotient
1074 C Jacobians.
1075 C
1076 C In the user-supplied version of DEWSET, it may be desirable to use
1077 C the current values of derivatives of y. Derivatives up to order NQ
1078 C are available from the history array YH, described above under
1079 C optional outputs. In DEWSET, YH is identical to the YCUR array,
1080 C extended to NQ + 1 columns with a column length of NYH and scale
1081 C factors of H**j/factorial(j). On the first call for the problem,
1082 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1083 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
1084 C can be obtained by including in SEWSET the statements:
1085 C DOUBLE PRECISION RLS
1086 C COMMON /DLS001/ RLS(218),ILS(37)
1087 C NQ = ILS(33)
1088 C NST = ILS(34)
1089 C H = RLS(212)
1090 C Thus, for example, the current value of dy/dt can be obtained as
1091 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary
1092 C when NST = 0).
1093 C
1094 C DVNORM
1095 C ------
1096 C DVNORM is a real function routine which computes the weighted
1097 C root-mean-square norm of a vector v:
1098 C
1099 C d = DVNORM (n, v, w)
1100 C
1101 C where:
1102 C n = the length of the vector,
1103 C v = real array of length n containing the vector,
1104 C w = real array of length n containing weights,
1105 C d = SQRT( (1/n) * sum(v(i)*w(i))**2 ).
1106 C
1107 C DVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where
1108 C EWT is as set by subroutine DEWSET.
1109 C
1110 C If the user supplies this function, it should return a nonnegative
1111 C value of DVNORM suitable for use in the error control in DLSODE.
1112 C None of the arguments should be altered by DVNORM. For example, a
1113 C user-supplied DVNORM routine might:
1114 C - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
1115 C - Ignore some components of v in the norm, with the effect of
1116 C suppressing the error control on those components of Y.
1117 C ---------------------------------------------------------------------
1118 C***ROUTINES CALLED DEWSET, DINTDY, DUMACH, DSTODE, DVNORM, XERRWD
1119 C***COMMON BLOCKS DLS001
1120 C***REVISION HISTORY (YYYYMMDD)
1121 C 19791129 DATE WRITTEN
1122 C 19791213 Minor changes to declarations; DELP init. in STODE.
1123 C 19800118 Treat NEQ as array; integer declarations added throughout;
1124 C minor changes to prologue.
1125 C 19800306 Corrected TESCO(1,NQP1) setting in CFODE.
1126 C 19800519 Corrected access of YH on forced order reduction;
1127 C numerous corrections to prologues and other comments.
1128 C 19800617 In main driver, added loading of SQRT(UROUND) in RWORK;
1129 C minor corrections to main prologue.
1130 C 19800923 Added zero initialization of HU and NQU.
1131 C 19801218 Revised XERRWD routine; minor corrections to main prologue.
1132 C 19810401 Minor changes to comments and an error message.
1133 C 19810814 Numerous revisions: replaced EWT by 1/EWT; used flags
1134 C JCUR, ICF, IERPJ, IERSL between STODE and subordinates;
1135 C added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF;
1136 C reorganized returns from STODE; reorganized type decls.;
1137 C fixed message length in XERRWD; changed default LUNIT to 6;
1138 C changed Common lengths; changed comments throughout.
1139 C 19870330 Major update by ACH: corrected comments throughout;
1140 C removed TRET from Common; rewrote EWSET with 4 loops;
1141 C fixed t test in INTDY; added Cray directives in STODE;
1142 C in STODE, fixed DELP init. and logic around PJAC call;
1143 C combined routines to save/restore Common;
1144 C passed LEVEL = 0 in error message calls (except run abort).
1145 C 19890426 Modified prologue to SLATEC/LDOC format. (FNF)
1146 C 19890501 Many improvements to prologue. (FNF)
1147 C 19890503 A few final corrections to prologue. (FNF)
1148 C 19890504 Minor cosmetic changes. (FNF)
1149 C 19890510 Corrected description of Y in Arguments section. (FNF)
1150 C 19890517 Minor corrections to prologue. (FNF)
1151 C 19920514 Updated with prologue edited 891025 by G. Shaw for manual.
1152 C 19920515 Converted source lines to upper case. (FNF)
1153 C 19920603 Revised XERRWD calls using mixed upper-lower case. (ACH)
1154 C 19920616 Revised prologue comment regarding CFT. (ACH)
1155 C 19921116 Revised prologue comments regarding Common. (ACH).
1156 C 19930326 Added comment about non-reentrancy. (FNF)
1157 C 19930723 Changed D1MACH to DUMACH. (FNF)
1158 C 19930801 Removed ILLIN and NTREP from Common (affects driver logic);
1159 C minor changes to prologue and internal comments;
1160 C changed Hollerith strings to quoted strings;
1161 C changed internal comments to mixed case;
1162 C replaced XERRWD with new version using character type;
1163 C changed dummy dimensions from 1 to *. (ACH)
1164 C 19930809 Changed to generic intrinsic names; changed names of
1165 C subprograms and Common blocks to DLSODE etc. (ACH)
1166 C 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH)
1167 C 20010412 Removed all 'own' variables from Common block /DLS001/
1168 C (affects declarations in 6 routines). (ACH)
1169 C 20010509 Minor corrections to prologue. (ACH)
1170 C 20031105 Restored 'own' variables to Common block /DLS001/, to
1171 C enable interrupt/restart feature. (ACH)
1172 C 20031112 Added SAVE statements for data-loaded constants.
1173 C
1174 C***END PROLOGUE DLSODE
1175 C
1176 C*Internal Notes:
1177 C
1178 C Other Routines in the DLSODE Package.
1179 C
1180 C In addition to Subroutine DLSODE, the DLSODE package includes the
1181 C following subroutines and function routines:
1182 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1183 C DSTODE is the core integrator, which does one step of the
1184 C integration and the associated error control.
1185 C DCFODE sets all method coefficients and test constants.
1186 C DPREPJ computes and preprocesses the Jacobian matrix J = df/dy
1187 C and the Newton iteration matrix P = I - h*l0*J.
1188 C DSOLSY manages solution of linear system in chord iteration.
1189 C DEWSET sets the error weight vector EWT before each step.
1190 C DVNORM computes the weighted R.M.S. norm of a vector.
1191 C DSRCOM is a user-callable routine to save and restore
1192 C the contents of the internal Common block.
1193 C DGEFA and DGESL are routines from LINPACK for solving full
1194 C systems of linear algebraic equations.
1195 C DGBFA and DGBSL are routines from LINPACK for solving banded
1196 C linear systems.
1197 C DUMACH computes the unit roundoff in a machine-independent manner.
1198 C XERRWD, XSETUN, XSETF, IXSAV, IUMACH handle the printing of all
1199 C error messages and warnings. XERRWD is machine-dependent.
1200 C Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1201 C All the others are subroutines.
1202 C
1203 C**End
1204 C
1205 C Declare externals.
1208 C
1209 C Declare all other variables.
1221 LOGICAL IHIT
1222 CHARACTER*80 MSG
1224 C-----------------------------------------------------------------------
1225 C The following internal Common block contains
1226 C (a) variables which are local to any subroutine but whose values must
1227 C be preserved between calls to the routine ("own" variables), and
1228 C (b) variables which are communicated between subroutines.
1229 C The block DLS001 is declared in subroutines DLSODE, DINTDY, DSTODE,
1230 C DPREPJ, and DSOLSY.
1231 C Groups of variables are replaced by dummy arrays in the Common
1232 C declarations in routines where those variables are not used.
1233 C-----------------------------------------------------------------------
1240 C
1242 C-----------------------------------------------------------------------
1243 C Block A.
1244 C This code block is executed on every call.
1245 C It tests ISTATE and ITASK for legality and branches appropriately.
1246 C If ISTATE .GT. 1 but the flag INIT shows that initialization has
1247 C not yet been done, an error return occurs.
1248 C If ISTATE = 1 and TOUT = T, return immediately.
1249 C-----------------------------------------------------------------------
1250 C
1251 C***FIRST EXECUTABLE STATEMENT DLSODE
1260 C-----------------------------------------------------------------------
1261 C Block B.
1262 C The next code block is executed for the initial call (ISTATE = 1),
1263 C or for a continuation call with parameter changes (ISTATE = 3).
1264 C It contains checking of all inputs and various initializations.
1265 C
1266 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1267 C MF, ML, and MU.
1268 C-----------------------------------------------------------------------
1285 C Next process and check the optional inputs. --------------------------
1288 MXSTEP = MXSTP0
1289 MXHNIL = MXHNL0
1313 C-----------------------------------------------------------------------
1314 C Set work array pointers and check lengths LRW and LIW.
1315 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1316 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1317 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1318 C-----------------------------------------------------------------------
1337 C Check RTOL and ATOL for legality. ------------------------------------
1347 C If ISTATE = 3, set flag to signal parameter changes to DSTODE. -------
1350 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1353 C Reload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
1356 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1363 C-----------------------------------------------------------------------
1364 C Block C.
1365 C The next block is for the initial call only (ISTATE = 1).
1366 C It contains all remaining initializations, the initial call to F,
1367 C and the calculation of the initial step size.
1368 C The error weights in EWT are inverted after being loaded.
1369 C-----------------------------------------------------------------------
1371 TN = T
1389 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1393 C Load the initial value vector in YH. ---------------------------------
1396 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1403 C-----------------------------------------------------------------------
1404 C The coding below computes the step size, H0, to be attempted on the
1405 C first step, unless the user has supplied a value for this.
1406 C First check that TOUT - T differs significantly from zero.
1407 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(I))
1408 C if this is positive, or MAX(ATOL(I)/ABS(Y(I))) otherwise, adjusted
1409 C so as to be between 100*UROUND and 1.0E-3.
1410 C Then the computed value H0 is given by..
1411 C NEQ
1412 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 )
1413 C 1
1414 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1415 C f(i) = i-th component of initial value of f,
1416 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1417 C The sign of H0 is inferred from the initial values of TOUT and T.
1418 C-----------------------------------------------------------------------
1441 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1444 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1445 H = H0
1449 C-----------------------------------------------------------------------
1450 C Block D.
1451 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1452 C and is to check stop conditions before taking a step.
1453 C-----------------------------------------------------------------------
1459 T = TOUT
1471 T = TOUT
1482 C-----------------------------------------------------------------------
1483 C Block E.
1484 C The next block is normally executed for all calls and contains
1485 C the call to the one-step core integrator DSTODE.
1486 C
1487 C This is a looping point for the integration steps.
1488 C
1489 C First check for too many steps being taken, update EWT (if not at
1490 C start of problem), check for too much accuracy being requested, and
1491 C check for H below the roundoff level in T.
1492 C-----------------------------------------------------------------------
1519 C-----------------------------------------------------------------------
1520 C CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPREPJ,DSOLSY)
1521 C-----------------------------------------------------------------------
1527 C-----------------------------------------------------------------------
1528 C Block F.
1529 C The following block handles the case of a successful return from the
1530 C core integrator (KFLAG = 0). Test for stop conditions.
1531 C-----------------------------------------------------------------------
1534 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1537 T = TOUT
1539 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1542 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1545 T = TOUT
1555 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1558 C-----------------------------------------------------------------------
1559 C Block G.
1560 C The following block handles all successful returns from DLSODE.
1561 C If ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
1562 C ISTATE is set to 2, and the optional outputs are loaded into the
1563 C work arrays before returning.
1564 C-----------------------------------------------------------------------
1567 T = TN
1579 RETURN
1580 C-----------------------------------------------------------------------
1581 C Block H.
1582 C The following block handles all unsuccessful returns other than
1583 C those for illegal input. First the error message routine is called.
1584 C If there was an error test or convergence test failure, IMXER is set.
1585 C Then Y is loaded from YH and T is set to TN. The optional outputs
1586 C are loaded into the work arrays before returning.
1587 C-----------------------------------------------------------------------
1588 C The maximum number of steps was taken before reaching TOUT. ----------
1595 C EWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
1601 C Too much accuracy requested for machine precision. -------------------
1609 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1616 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1624 C Compute IMXER if relevant. -------------------------------------------
1630 BIG = SIZE
1631 IMXER = I
1634 C Set Y vector, T, and optional outputs. -------------------------------
1637 T = TN
1646 RETURN
1647 C-----------------------------------------------------------------------
1648 C Block I.
1649 C The following block handles all error returns due to illegal input
1650 C (ISTATE = -3), as detected before calling the core integrator.
1651 C First the error message routine is called. If the illegal input
1652 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1653 C-----------------------------------------------------------------------
1747 C
1749 RETURN
1750 C
1753 RETURN
1754 C----------------------- END OF SUBROUTINE DLSODE ----------------------
1755 END