| blob | f1adf2f4f81ffb93ad89b25fae0db80f582d7863 |
1 *DECK DLSODA
8 C-----------------------------------------------------------------------
9 C This is the 12 November 2003 version of
10 C DLSODA: Livermore Solver for Ordinary Differential Equations, with
11 C Automatic method switching for stiff and nonstiff problems.
12 C
13 C This version is in double precision.
14 C
15 C DLSODA solves the initial value problem for stiff or nonstiff
16 C systems of first order ODEs,
17 C dy/dt = f(t,y) , or, in component form,
18 C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
19 C
20 C This a variant version of the DLSODE package.
21 C It switches automatically between stiff and nonstiff methods.
22 C This means that the user does not have to determine whether the
23 C problem is stiff or not, and the solver will automatically choose the
24 C appropriate method. It always starts with the nonstiff method.
25 C
26 C Authors: Alan C. Hindmarsh
27 C Center for Applied Scientific Computing, L-561
28 C Lawrence Livermore National Laboratory
29 C Livermore, CA 94551
30 C and
31 C Linda R. Petzold
32 C Univ. of California at Santa Barbara
33 C Dept. of Computer Science
34 C Santa Barbara, CA 93106
35 C
36 C References:
37 C 1. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE
38 C Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.),
39 C North-Holland, Amsterdam, 1983, pp. 55-64.
40 C 2. Linda R. Petzold, Automatic Selection of Methods for Solving
41 C Stiff and Nonstiff Systems of Ordinary Differential Equations,
42 C Siam J. Sci. Stat. Comput. 4 (1983), pp. 136-148.
43 C-----------------------------------------------------------------------
44 C Summary of Usage.
45 C
46 C Communication between the user and the DLSODA package, for normal
47 C situations, is summarized here. This summary describes only a subset
48 C of the full set of options available. See the full description for
49 C details, including alternative treatment of the Jacobian matrix,
50 C optional inputs and outputs, nonstandard options, and
51 C instructions for special situations. See also the example
52 C problem (with program and output) following this summary.
53 C
54 C A. First provide a subroutine of the form:
55 C SUBROUTINE F (NEQ, T, Y, YDOT)
56 C DOUBLE PRECISION T, Y(*), YDOT(*)
57 C which supplies the vector function f by loading YDOT(i) with f(i).
58 C
59 C B. Write a main program which calls Subroutine DLSODA once for
60 C each point at which answers are desired. This should also provide
61 C for possible use of logical unit 6 for output of error messages
62 C by DLSODA. On the first call to DLSODA, supply arguments as follows:
63 C F = name of subroutine for right-hand side vector f.
64 C This name must be declared External in calling program.
65 C NEQ = number of first order ODEs.
66 C Y = array of initial values, of length NEQ.
67 C T = the initial value of the independent variable.
68 C TOUT = first point where output is desired (.ne. T).
69 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
70 C RTOL = relative tolerance parameter (scalar).
71 C ATOL = absolute tolerance parameter (scalar or array).
72 C the estimated local error in y(i) will be controlled so as
73 C to be less than
74 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
75 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
76 C Thus the local error test passes if, in each component,
77 C either the absolute error is less than ATOL (or ATOL(i)),
78 C or the relative error is less than RTOL.
79 C Use RTOL = 0.0 for pure absolute error control, and
80 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
81 C control. Caution: actual (global) errors may exceed these
82 C local tolerances, so choose them conservatively.
83 C ITASK = 1 for normal computation of output values of y at t = TOUT.
84 C ISTATE = integer flag (input and output). Set ISTATE = 1.
85 C IOPT = 0 to indicate no optional inputs used.
86 C RWORK = real work array of length at least:
87 C 22 + NEQ * MAX(16, NEQ + 9).
88 C See also Paragraph E below.
89 C LRW = declared length of RWORK (in user's dimension).
90 C IWORK = integer work array of length at least 20 + NEQ.
91 C LIW = declared length of IWORK (in user's dimension).
92 C JAC = name of subroutine for Jacobian matrix.
93 C Use a dummy name. See also Paragraph E below.
94 C JT = Jacobian type indicator. Set JT = 2.
95 C See also Paragraph E below.
96 C Note that the main program must declare arrays Y, RWORK, IWORK,
97 C and possibly ATOL.
98 C
99 C C. The output from the first call (or any call) is:
100 C Y = array of computed values of y(t) vector.
101 C T = corresponding value of independent variable (normally TOUT).
102 C ISTATE = 2 if DLSODA was successful, negative otherwise.
103 C -1 means excess work done on this call (perhaps wrong JT).
104 C -2 means excess accuracy requested (tolerances too small).
105 C -3 means illegal input detected (see printed message).
106 C -4 means repeated error test failures (check all inputs).
107 C -5 means repeated convergence failures (perhaps bad Jacobian
108 C supplied or wrong choice of JT or tolerances).
109 C -6 means error weight became zero during problem. (Solution
110 C component i vanished, and ATOL or ATOL(i) = 0.)
111 C -7 means work space insufficient to finish (see messages).
112 C
113 C D. To continue the integration after a successful return, simply
114 C reset TOUT and call DLSODA again. No other parameters need be reset.
115 C
116 C E. Note: If and when DLSODA regards the problem as stiff, and
117 C switches methods accordingly, it must make use of the NEQ by NEQ
118 C Jacobian matrix, J = df/dy. For the sake of simplicity, the
119 C inputs to DLSODA recommended in Paragraph B above cause DLSODA to
120 C treat J as a full matrix, and to approximate it internally by
121 C difference quotients. Alternatively, J can be treated as a band
122 C matrix (with great potential reduction in the size of the RWORK
123 C array). Also, in either the full or banded case, the user can supply
124 C J in closed form, with a routine whose name is passed as the JAC
125 C argument. These alternatives are described in the paragraphs on
126 C RWORK, JAC, and JT in the full description of the call sequence below.
127 C
128 C-----------------------------------------------------------------------
129 C Example Problem.
130 C
131 C The following is a simple example problem, with the coding
132 C needed for its solution by DLSODA. The problem is from chemical
133 C kinetics, and consists of the following three rate equations:
134 C dy1/dt = -.04*y1 + 1.e4*y2*y3
135 C dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
136 C dy3/dt = 3.e7*y2**2
137 C on the interval from t = 0.0 to t = 4.e10, with initial conditions
138 C y1 = 1.0, y2 = y3 = 0. The problem is stiff.
139 C
140 C The following coding solves this problem with DLSODA,
141 C printing results at t = .4, 4., ..., 4.e10. It uses
142 C ITOL = 2 and ATOL much smaller for y2 than y1 or y3 because
143 C y2 has much smaller values.
144 C At the end of the run, statistical quantities of interest are
145 C printed (see optional outputs in the full description below).
146 C
147 C EXTERNAL FEX
148 C DOUBLE PRECISION ATOL, RTOL, RWORK, T, TOUT, Y
149 C DIMENSION Y(3), ATOL(3), RWORK(70), IWORK(23)
150 C NEQ = 3
151 C Y(1) = 1.
152 C Y(2) = 0.
153 C Y(3) = 0.
154 C T = 0.
155 C TOUT = .4
156 C ITOL = 2
157 C RTOL = 1.D-4
158 C ATOL(1) = 1.D-6
159 C ATOL(2) = 1.D-10
160 C ATOL(3) = 1.D-6
161 C ITASK = 1
162 C ISTATE = 1
163 C IOPT = 0
164 C LRW = 70
165 C LIW = 23
166 C JT = 2
167 C DO 40 IOUT = 1,12
168 C CALL DLSODA(FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,
169 C 1 IOPT,RWORK,LRW,IWORK,LIW,JDUM,JT)
170 C WRITE(6,20)T,Y(1),Y(2),Y(3)
171 C 20 FORMAT(' At t =',D12.4,' Y =',3D14.6)
172 C IF (ISTATE .LT. 0) GO TO 80
173 C 40 TOUT = TOUT*10.
174 C WRITE(6,60)IWORK(11),IWORK(12),IWORK(13),IWORK(19),RWORK(15)
175 C 60 FORMAT(/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4/
176 C 1 ' Method last used =',I2,' Last switch was at t =',D12.4)
177 C STOP
178 C 80 WRITE(6,90)ISTATE
179 C 90 FORMAT(///' Error halt.. ISTATE =',I3)
180 C STOP
181 C END
182 C
183 C SUBROUTINE FEX (NEQ, T, Y, YDOT)
184 C DOUBLE PRECISION T, Y, YDOT
185 C DIMENSION Y(3), YDOT(3)
186 C YDOT(1) = -.04*Y(1) + 1.D4*Y(2)*Y(3)
187 C YDOT(3) = 3.D7*Y(2)*Y(2)
188 C YDOT(2) = -YDOT(1) - YDOT(3)
189 C RETURN
190 C END
191 C
192 C The output of this program (on a CDC-7600 in single precision)
193 C is as follows:
194 C
195 C At t = 4.0000e-01 y = 9.851712e-01 3.386380e-05 1.479493e-02
196 C At t = 4.0000e+00 Y = 9.055333e-01 2.240655e-05 9.444430e-02
197 C At t = 4.0000e+01 Y = 7.158403e-01 9.186334e-06 2.841505e-01
198 C At t = 4.0000e+02 Y = 4.505250e-01 3.222964e-06 5.494717e-01
199 C At t = 4.0000e+03 Y = 1.831975e-01 8.941774e-07 8.168016e-01
200 C At t = 4.0000e+04 Y = 3.898730e-02 1.621940e-07 9.610125e-01
201 C At t = 4.0000e+05 Y = 4.936363e-03 1.984221e-08 9.950636e-01
202 C At t = 4.0000e+06 Y = 5.161831e-04 2.065786e-09 9.994838e-01
203 C At t = 4.0000e+07 Y = 5.179817e-05 2.072032e-10 9.999482e-01
204 C At t = 4.0000e+08 Y = 5.283401e-06 2.113371e-11 9.999947e-01
205 C At t = 4.0000e+09 Y = 4.659031e-07 1.863613e-12 9.999995e-01
206 C At t = 4.0000e+10 Y = 1.404280e-08 5.617126e-14 1.000000e+00
207 C
208 C No. steps = 361 No. f-s = 693 No. J-s = 64
209 C Method last used = 2 Last switch was at t = 6.0092e-03
210 C-----------------------------------------------------------------------
211 C Full description of user interface to DLSODA.
212 C
213 C The user interface to DLSODA consists of the following parts.
214 C
215 C 1. The call sequence to Subroutine DLSODA, which is a driver
216 C routine for the solver. This includes descriptions of both
217 C the call sequence arguments and of user-supplied routines.
218 C following these descriptions is a description of
219 C optional inputs available through the call sequence, and then
220 C a description of optional outputs (in the work arrays).
221 C
222 C 2. Descriptions of other routines in the DLSODA package that may be
223 C (optionally) called by the user. These provide the ability to
224 C alter error message handling, save and restore the internal
225 C Common, and obtain specified derivatives of the solution y(t).
226 C
227 C 3. Descriptions of Common blocks to be declared in overlay
228 C or similar environments, or to be saved when doing an interrupt
229 C of the problem and continued solution later.
230 C
231 C 4. Description of a subroutine in the DLSODA package,
232 C which the user may replace with his/her own version, if desired.
233 C this relates to the measurement of errors.
234 C
235 C-----------------------------------------------------------------------
236 C Part 1. Call Sequence.
237 C
238 C The call sequence parameters used for input only are
239 C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, JT,
240 C and those used for both input and output are
241 C Y, T, ISTATE.
242 C The work arrays RWORK and IWORK are also used for conditional and
243 C optional inputs and optional outputs. (The term output here refers
244 C to the return from Subroutine DLSODA to the user's calling program.)
245 C
246 C The legality of input parameters will be thoroughly checked on the
247 C initial call for the problem, but not checked thereafter unless a
248 C change in input parameters is flagged by ISTATE = 3 on input.
249 C
250 C The descriptions of the call arguments are as follows.
251 C
252 C F = the name of the user-supplied subroutine defining the
253 C ODE system. The system must be put in the first-order
254 C form dy/dt = f(t,y), where f is a vector-valued function
255 C of the scalar t and the vector y. Subroutine F is to
256 C compute the function f. It is to have the form
257 C SUBROUTINE F (NEQ, T, Y, YDOT)
258 C DOUBLE PRECISION T, Y(*), YDOT(*)
259 C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
260 C is output. Y and YDOT are arrays of length NEQ.
261 C Subroutine F should not alter Y(1),...,Y(NEQ).
262 C F must be declared External in the calling program.
263 C
264 C Subroutine F may access user-defined quantities in
265 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
266 C (dimensioned in F) and/or Y has length exceeding NEQ(1).
267 C See the descriptions of NEQ and Y below.
268 C
269 C If quantities computed in the F routine are needed
270 C externally to DLSODA, an extra call to F should be made
271 C for this purpose, for consistent and accurate results.
272 C If only the derivative dy/dt is needed, use DINTDY instead.
273 C
274 C NEQ = the size of the ODE system (number of first order
275 C ordinary differential equations). Used only for input.
276 C NEQ may be decreased, but not increased, during the problem.
277 C If NEQ is decreased (with ISTATE = 3 on input), the
278 C remaining components of Y should be left undisturbed, if
279 C these are to be accessed in F and/or JAC.
280 C
281 C Normally, NEQ is a scalar, and it is generally referred to
282 C as a scalar in this user interface description. However,
283 C NEQ may be an array, with NEQ(1) set to the system size.
284 C (The DLSODA package accesses only NEQ(1).) In either case,
285 C this parameter is passed as the NEQ argument in all calls
286 C to F and JAC. Hence, if it is an array, locations
287 C NEQ(2),... may be used to store other integer data and pass
288 C it to F and/or JAC. Subroutines F and/or JAC must include
289 C NEQ in a Dimension statement in that case.
290 C
291 C Y = a real array for the vector of dependent variables, of
292 C length NEQ or more. Used for both input and output on the
293 C first call (ISTATE = 1), and only for output on other calls.
294 C On the first call, Y must contain the vector of initial
295 C values. On output, Y contains the computed solution vector,
296 C evaluated at T. If desired, the Y array may be used
297 C for other purposes between calls to the solver.
298 C
299 C This array is passed as the Y argument in all calls to
300 C F and JAC. Hence its length may exceed NEQ, and locations
301 C Y(NEQ+1),... may be used to store other real data and
302 C pass it to F and/or JAC. (The DLSODA package accesses only
303 C Y(1),...,Y(NEQ).)
304 C
305 C T = the independent variable. On input, T is used only on the
306 C first call, as the initial point of the integration.
307 C on output, after each call, T is the value at which a
308 C computed solution Y is evaluated (usually the same as TOUT).
309 C on an error return, T is the farthest point reached.
310 C
311 C TOUT = the next value of t at which a computed solution is desired.
312 C Used only for input.
313 C
314 C When starting the problem (ISTATE = 1), TOUT may be equal
315 C to T for one call, then should .ne. T for the next call.
316 C For the initial t, an input value of TOUT .ne. T is used
317 C in order to determine the direction of the integration
318 C (i.e. the algebraic sign of the step sizes) and the rough
319 C scale of the problem. Integration in either direction
320 C (forward or backward in t) is permitted.
321 C
322 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
323 C the first call (i.e. the first call with TOUT .ne. T).
324 C Otherwise, TOUT is required on every call.
325 C
326 C If ITASK = 1, 3, or 4, the values of TOUT need not be
327 C monotone, but a value of TOUT which backs up is limited
328 C to the current internal T interval, whose endpoints are
329 C TCUR - HU and TCUR (see optional outputs, below, for
330 C TCUR and HU).
331 C
332 C ITOL = an indicator for the type of error control. See
333 C description below under ATOL. Used only for input.
334 C
335 C RTOL = a relative error tolerance parameter, either a scalar or
336 C an array of length NEQ. See description below under ATOL.
337 C Input only.
338 C
339 C ATOL = an absolute error tolerance parameter, either a scalar or
340 C an array of length NEQ. Input only.
341 C
342 C The input parameters ITOL, RTOL, and ATOL determine
343 C the error control performed by the solver. The solver will
344 C control the vector E = (E(i)) of estimated local errors
345 C in y, according to an inequality of the form
346 C max-norm of ( E(i)/EWT(i) ) .le. 1,
347 C where EWT = (EWT(i)) is a vector of positive error weights.
348 C The values of RTOL and ATOL should all be non-negative.
349 C The following table gives the types (scalar/array) of
350 C RTOL and ATOL, and the corresponding form of EWT(i).
351 C
352 C ITOL RTOL ATOL EWT(i)
353 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
354 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
355 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
356 C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i)
357 C
358 C When either of these parameters is a scalar, it need not
359 C be dimensioned in the user's calling program.
360 C
361 C If none of the above choices (with ITOL, RTOL, and ATOL
362 C fixed throughout the problem) is suitable, more general
363 C error controls can be obtained by substituting a
364 C user-supplied routine for the setting of EWT.
365 C See Part 4 below.
366 C
367 C If global errors are to be estimated by making a repeated
368 C run on the same problem with smaller tolerances, then all
369 C components of RTOL and ATOL (i.e. of EWT) should be scaled
370 C down uniformly.
371 C
372 C ITASK = an index specifying the task to be performed.
373 C Input only. ITASK has the following values and meanings.
374 C 1 means normal computation of output values of y(t) at
375 C t = TOUT (by overshooting and interpolating).
376 C 2 means take one step only and return.
377 C 3 means stop at the first internal mesh point at or
378 C beyond t = TOUT and return.
379 C 4 means normal computation of output values of y(t) at
380 C t = TOUT but without overshooting t = TCRIT.
381 C TCRIT must be input as RWORK(1). TCRIT may be equal to
382 C or beyond TOUT, but not behind it in the direction of
383 C integration. This option is useful if the problem
384 C has a singularity at or beyond t = TCRIT.
385 C 5 means take one step, without passing TCRIT, and return.
386 C TCRIT must be input as RWORK(1).
387 C
388 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
389 C (within roundoff), it will return T = TCRIT (exactly) to
390 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
391 C in which case answers at t = TOUT are returned first).
392 C
393 C ISTATE = an index used for input and output to specify the
394 C the state of the calculation.
395 C
396 C On input, the values of ISTATE are as follows.
397 C 1 means this is the first call for the problem
398 C (initializations will be done). See note below.
399 C 2 means this is not the first call, and the calculation
400 C is to continue normally, with no change in any input
401 C parameters except possibly TOUT and ITASK.
402 C (If ITOL, RTOL, and/or ATOL are changed between calls
403 C with ISTATE = 2, the new values will be used but not
404 C tested for legality.)
405 C 3 means this is not the first call, and the
406 C calculation is to continue normally, but with
407 C a change in input parameters other than
408 C TOUT and ITASK. Changes are allowed in
409 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, JT, ML, MU,
410 C and any optional inputs except H0, MXORDN, and MXORDS.
411 C (See IWORK description for ML and MU.)
412 C Note: A preliminary call with TOUT = T is not counted
413 C as a first call here, as no initialization or checking of
414 C input is done. (Such a call is sometimes useful for the
415 C purpose of outputting the initial conditions.)
416 C Thus the first call for which TOUT .ne. T requires
417 C ISTATE = 1 on input.
418 C
419 C On output, ISTATE has the following values and meanings.
420 C 1 means nothing was done; TOUT = T and ISTATE = 1 on input.
421 C 2 means the integration was performed successfully.
422 C -1 means an excessive amount of work (more than MXSTEP
423 C steps) was done on this call, before completing the
424 C requested task, but the integration was otherwise
425 C successful as far as T. (MXSTEP is an optional input
426 C and is normally 500.) To continue, the user may
427 C simply reset ISTATE to a value .gt. 1 and call again
428 C (the excess work step counter will be reset to 0).
429 C In addition, the user may increase MXSTEP to avoid
430 C this error return (see below on optional inputs).
431 C -2 means too much accuracy was requested for the precision
432 C of the machine being used. This was detected before
433 C completing the requested task, but the integration
434 C was successful as far as T. To continue, the tolerance
435 C parameters must be reset, and ISTATE must be set
436 C to 3. The optional output TOLSF may be used for this
437 C purpose. (Note: If this condition is detected before
438 C taking any steps, then an illegal input return
439 C (ISTATE = -3) occurs instead.)
440 C -3 means illegal input was detected, before taking any
441 C integration steps. See written message for details.
442 C Note: If the solver detects an infinite loop of calls
443 C to the solver with illegal input, it will cause
444 C the run to stop.
445 C -4 means there were repeated error test failures on
446 C one attempted step, before completing the requested
447 C task, but the integration was successful as far as T.
448 C The problem may have a singularity, or the input
449 C may be inappropriate.
450 C -5 means there were repeated convergence test failures on
451 C one attempted step, before completing the requested
452 C task, but the integration was successful as far as T.
453 C This may be caused by an inaccurate Jacobian matrix,
454 C if one is being used.
455 C -6 means EWT(i) became zero for some i during the
456 C integration. Pure relative error control (ATOL(i)=0.0)
457 C was requested on a variable which has now vanished.
458 C The integration was successful as far as T.
459 C -7 means the length of RWORK and/or IWORK was too small to
460 C proceed, but the integration was successful as far as T.
461 C This happens when DLSODA chooses to switch methods
462 C but LRW and/or LIW is too small for the new method.
463 C
464 C Note: Since the normal output value of ISTATE is 2,
465 C it does not need to be reset for normal continuation.
466 C Also, since a negative input value of ISTATE will be
467 C regarded as illegal, a negative output value requires the
468 C user to change it, and possibly other inputs, before
469 C calling the solver again.
470 C
471 C IOPT = an integer flag to specify whether or not any optional
472 C inputs are being used on this call. Input only.
473 C The optional inputs are listed separately below.
474 C IOPT = 0 means no optional inputs are being used.
475 C default values will be used in all cases.
476 C IOPT = 1 means one or more optional inputs are being used.
477 C
478 C RWORK = a real array (double precision) for work space, and (in the
479 C first 20 words) for conditional and optional inputs and
480 C optional outputs.
481 C As DLSODA switches automatically between stiff and nonstiff
482 C methods, the required length of RWORK can change during the
483 C problem. Thus the RWORK array passed to DLSODA can either
484 C have a static (fixed) length large enough for both methods,
485 C or have a dynamic (changing) length altered by the calling
486 C program in response to output from DLSODA.
487 C
488 C --- Fixed Length Case ---
489 C If the RWORK length is to be fixed, it should be at least
490 C MAX (LRN, LRS),
491 C where LRN and LRS are the RWORK lengths required when the
492 C current method is nonstiff or stiff, respectively.
493 C
494 C The separate RWORK length requirements LRN and LRS are
495 C as follows:
496 C IF NEQ is constant and the maximum method orders have
497 C their default values, then
498 C LRN = 20 + 16*NEQ,
499 C LRS = 22 + 9*NEQ + NEQ**2 if JT = 1 or 2,
500 C LRS = 22 + 10*NEQ + (2*ML+MU)*NEQ if JT = 4 or 5.
501 C Under any other conditions, LRN and LRS are given by:
502 C LRN = 20 + NYH*(MXORDN+1) + 3*NEQ,
503 C LRS = 20 + NYH*(MXORDS+1) + 3*NEQ + LMAT,
504 C where
505 C NYH = the initial value of NEQ,
506 C MXORDN = 12, unless a smaller value is given as an
507 C optional input,
508 C MXORDS = 5, unless a smaller value is given as an
509 C optional input,
510 C LMAT = length of matrix work space:
511 C LMAT = NEQ**2 + 2 if JT = 1 or 2,
512 C LMAT = (2*ML + MU + 1)*NEQ + 2 if JT = 4 or 5.
513 C
514 C --- Dynamic Length Case ---
515 C If the length of RWORK is to be dynamic, then it should
516 C be at least LRN or LRS, as defined above, depending on the
517 C current method. Initially, it must be at least LRN (since
518 C DLSODA starts with the nonstiff method). On any return
519 C from DLSODA, the optional output MCUR indicates the current
520 C method. If MCUR differs from the value it had on the
521 C previous return, or if there has only been one call to
522 C DLSODA and MCUR is now 2, then DLSODA has switched
523 C methods during the last call, and the length of RWORK
524 C should be reset (to LRN if MCUR = 1, or to LRS if
525 C MCUR = 2). (An increase in the RWORK length is required
526 C if DLSODA returned ISTATE = -7, but not otherwise.)
527 C After resetting the length, call DLSODA with ISTATE = 3
528 C to signal that change.
529 C
530 C LRW = the length of the array RWORK, as declared by the user.
531 C (This will be checked by the solver.)
532 C
533 C IWORK = an integer array for work space.
534 C As DLSODA switches automatically between stiff and nonstiff
535 C methods, the required length of IWORK can change during
536 C problem, between
537 C LIS = 20 + NEQ and LIN = 20,
538 C respectively. Thus the IWORK array passed to DLSODA can
539 C either have a fixed length of at least 20 + NEQ, or have a
540 C dynamic length of at least LIN or LIS, depending on the
541 C current method. The comments on dynamic length under
542 C RWORK above apply here. Initially, this length need
543 C only be at least LIN = 20.
544 C
545 C The first few words of IWORK are used for conditional and
546 C optional inputs and optional outputs.
547 C
548 C The following 2 words in IWORK are conditional inputs:
549 C IWORK(1) = ML these are the lower and upper
550 C IWORK(2) = MU half-bandwidths, respectively, of the
551 C banded Jacobian, excluding the main diagonal.
552 C The band is defined by the matrix locations
553 C (i,j) with i-ML .le. j .le. i+MU. ML and MU
554 C must satisfy 0 .le. ML,MU .le. NEQ-1.
555 C These are required if JT is 4 or 5, and
556 C ignored otherwise. ML and MU may in fact be
557 C the band parameters for a matrix to which
558 C df/dy is only approximately equal.
559 C
560 C LIW = the length of the array IWORK, as declared by the user.
561 C (This will be checked by the solver.)
562 C
563 C Note: The base addresses of the work arrays must not be
564 C altered between calls to DLSODA for the same problem.
565 C The contents of the work arrays must not be altered
566 C between calls, except possibly for the conditional and
567 C optional inputs, and except for the last 3*NEQ words of RWORK.
568 C The latter space is used for internal scratch space, and so is
569 C available for use by the user outside DLSODA between calls, if
570 C desired (but not for use by F or JAC).
571 C
572 C JAC = the name of the user-supplied routine to compute the
573 C Jacobian matrix, df/dy, if JT = 1 or 4. The JAC routine
574 C is optional, but if the problem is expected to be stiff much
575 C of the time, you are encouraged to supply JAC, for the sake
576 C of efficiency. (Alternatively, set JT = 2 or 5 to have
577 C DLSODA compute df/dy internally by difference quotients.)
578 C If and when DLSODA uses df/dy, it treats this NEQ by NEQ
579 C matrix either as full (JT = 1 or 2), or as banded (JT =
580 C 4 or 5) with half-bandwidths ML and MU (discussed under
581 C IWORK above). In either case, if JT = 1 or 4, the JAC
582 C routine must compute df/dy as a function of the scalar t
583 C and the vector y. It is to have the form
584 C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
585 C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
586 C where NEQ, T, Y, ML, MU, and NROWPD are input and the array
587 C PD is to be loaded with partial derivatives (elements of
588 C the Jacobian matrix) on output. PD must be given a first
589 C dimension of NROWPD. T and Y have the same meaning as in
590 C Subroutine F.
591 C In the full matrix case (JT = 1), ML and MU are
592 C ignored, and the Jacobian is to be loaded into PD in
593 C columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
594 C In the band matrix case (JT = 4), the elements
595 C within the band are to be loaded into PD in columnwise
596 C manner, with diagonal lines of df/dy loaded into the rows
597 C of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j).
598 C ML and MU are the half-bandwidth parameters (see IWORK).
599 C The locations in PD in the two triangular areas which
600 C correspond to nonexistent matrix elements can be ignored
601 C or loaded arbitrarily, as they are overwritten by DLSODA.
602 C JAC need not provide df/dy exactly. A crude
603 C approximation (possibly with a smaller bandwidth) will do.
604 C In either case, PD is preset to zero by the solver,
605 C so that only the nonzero elements need be loaded by JAC.
606 C Each call to JAC is preceded by a call to F with the same
607 C arguments NEQ, T, and Y. Thus to gain some efficiency,
608 C intermediate quantities shared by both calculations may be
609 C saved in a user Common block by F and not recomputed by JAC,
610 C if desired. Also, JAC may alter the Y array, if desired.
611 C JAC must be declared External in the calling program.
612 C Subroutine JAC may access user-defined quantities in
613 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
614 C (dimensioned in JAC) and/or Y has length exceeding NEQ(1).
615 C See the descriptions of NEQ and Y above.
616 C
617 C JT = Jacobian type indicator. Used only for input.
618 C JT specifies how the Jacobian matrix df/dy will be
619 C treated, if and when DLSODA requires this matrix.
620 C JT has the following values and meanings:
621 C 1 means a user-supplied full (NEQ by NEQ) Jacobian.
622 C 2 means an internally generated (difference quotient) full
623 C Jacobian (using NEQ extra calls to F per df/dy value).
624 C 4 means a user-supplied banded Jacobian.
625 C 5 means an internally generated banded Jacobian (using
626 C ML+MU+1 extra calls to F per df/dy evaluation).
627 C If JT = 1 or 4, the user must supply a Subroutine JAC
628 C (the name is arbitrary) as described above under JAC.
629 C If JT = 2 or 5, a dummy argument can be used.
630 C-----------------------------------------------------------------------
631 C Optional Inputs.
632 C
633 C The following is a list of the optional inputs provided for in the
634 C call sequence. (See also Part 2.) For each such input variable,
635 C this table lists its name as used in this documentation, its
636 C location in the call sequence, its meaning, and the default value.
637 C The use of any of these inputs requires IOPT = 1, and in that
638 C case all of these inputs are examined. A value of zero for any
639 C of these optional inputs will cause the default value to be used.
640 C Thus to use a subset of the optional inputs, simply preload
641 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
642 C then set those of interest to nonzero values.
643 C
644 C Name Location Meaning and Default Value
645 C
646 C H0 RWORK(5) the step size to be attempted on the first step.
647 C The default value is determined by the solver.
648 C
649 C HMAX RWORK(6) the maximum absolute step size allowed.
650 C The default value is infinite.
651 C
652 C HMIN RWORK(7) the minimum absolute step size allowed.
653 C The default value is 0. (This lower bound is not
654 C enforced on the final step before reaching TCRIT
655 C when ITASK = 4 or 5.)
656 C
657 C IXPR IWORK(5) flag to generate extra printing at method switches.
658 C IXPR = 0 means no extra printing (the default).
659 C IXPR = 1 means print data on each switch.
660 C T, H, and NST will be printed on the same logical
661 C unit as used for error messages.
662 C
663 C MXSTEP IWORK(6) maximum number of (internally defined) steps
664 C allowed during one call to the solver.
665 C The default value is 500.
666 C
667 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
668 C warning that T + H = T on a step (H = step size).
669 C This must be positive to result in a non-default
670 C value. The default value is 10.
671 C
672 C MXORDN IWORK(8) the maximum order to be allowed for the nonstiff
673 C (Adams) method. the default value is 12.
674 C if MXORDN exceeds the default value, it will
675 C be reduced to the default value.
676 C MXORDN is held constant during the problem.
677 C
678 C MXORDS IWORK(9) the maximum order to be allowed for the stiff
679 C (BDF) method. The default value is 5.
680 C If MXORDS exceeds the default value, it will
681 C be reduced to the default value.
682 C MXORDS is held constant during the problem.
683 C-----------------------------------------------------------------------
684 C Optional Outputs.
685 C
686 C As optional additional output from DLSODA, the variables listed
687 C below are quantities related to the performance of DLSODA
688 C which are available to the user. These are communicated by way of
689 C the work arrays, but also have internal mnemonic names as shown.
690 C except where stated otherwise, all of these outputs are defined
691 C on any successful return from DLSODA, and on any return with
692 C ISTATE = -1, -2, -4, -5, or -6. On an illegal input return
693 C (ISTATE = -3), they will be unchanged from their existing values
694 C (if any), except possibly for TOLSF, LENRW, and LENIW.
695 C On any error return, outputs relevant to the error will be defined,
696 C as noted below.
697 C
698 C Name Location Meaning
699 C
700 C HU RWORK(11) the step size in t last used (successfully).
701 C
702 C HCUR RWORK(12) the step size to be attempted on the next step.
703 C
704 C TCUR RWORK(13) the current value of the independent variable
705 C which the solver has actually reached, i.e. the
706 C current internal mesh point in t. On output, TCUR
707 C will always be at least as far as the argument
708 C T, but may be farther (if interpolation was done).
709 C
710 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
711 C computed when a request for too much accuracy was
712 C detected (ISTATE = -3 if detected at the start of
713 C the problem, ISTATE = -2 otherwise). If ITOL is
714 C left unaltered but RTOL and ATOL are uniformly
715 C scaled up by a factor of TOLSF for the next call,
716 C then the solver is deemed likely to succeed.
717 C (The user may also ignore TOLSF and alter the
718 C tolerance parameters in any other way appropriate.)
719 C
720 C TSW RWORK(15) the value of t at the time of the last method
721 C switch, if any.
722 C
723 C NST IWORK(11) the number of steps taken for the problem so far.
724 C
725 C NFE IWORK(12) the number of f evaluations for the problem so far.
726 C
727 C NJE IWORK(13) the number of Jacobian evaluations (and of matrix
728 C LU decompositions) for the problem so far.
729 C
730 C NQU IWORK(14) the method order last used (successfully).
731 C
732 C NQCUR IWORK(15) the order to be attempted on the next step.
733 C
734 C IMXER IWORK(16) the index of the component of largest magnitude in
735 C the weighted local error vector ( E(i)/EWT(i) ),
736 C on an error return with ISTATE = -4 or -5.
737 C
738 C LENRW IWORK(17) the length of RWORK actually required, assuming
739 C that the length of RWORK is to be fixed for the
740 C rest of the problem, and that switching may occur.
741 C This is defined on normal returns and on an illegal
742 C input return for insufficient storage.
743 C
744 C LENIW IWORK(18) the length of IWORK actually required, assuming
745 C that the length of IWORK is to be fixed for the
746 C rest of the problem, and that switching may occur.
747 C This is defined on normal returns and on an illegal
748 C input return for insufficient storage.
749 C
750 C MUSED IWORK(19) the method indicator for the last successful step:
751 C 1 means Adams (nonstiff), 2 means BDF (stiff).
752 C
753 C MCUR IWORK(20) the current method indicator:
754 C 1 means Adams (nonstiff), 2 means BDF (stiff).
755 C This is the method to be attempted
756 C on the next step. Thus it differs from MUSED
757 C only if a method switch has just been made.
758 C
759 C The following two arrays are segments of the RWORK array which
760 C may also be of interest to the user as optional outputs.
761 C For each array, the table below gives its internal name,
762 C its base address in RWORK, and its description.
763 C
764 C Name Base Address Description
765 C
766 C YH 21 the Nordsieck history array, of size NYH by
767 C (NQCUR + 1), where NYH is the initial value
768 C of NEQ. For j = 0,1,...,NQCUR, column j+1
769 C of YH contains HCUR**j/factorial(j) times
770 C the j-th derivative of the interpolating
771 C polynomial currently representing the solution,
772 C evaluated at T = TCUR.
773 C
774 C ACOR LACOR array of size NEQ used for the accumulated
775 C (from Common corrections on each step, scaled on output
776 C as noted) to represent the estimated local error in y
777 C on the last step. This is the vector E in
778 C the description of the error control. It is
779 C defined only on a successful return from
780 C DLSODA. The base address LACOR is obtained by
781 C including in the user's program the
782 C following 2 lines:
783 C COMMON /DLS001/ RLS(218), ILS(37)
784 C LACOR = ILS(22)
785 C
786 C-----------------------------------------------------------------------
787 C Part 2. Other Routines Callable.
788 C
789 C The following are optional calls which the user may make to
790 C gain additional capabilities in conjunction with DLSODA.
791 C (The routines XSETUN and XSETF are designed to conform to the
792 C SLATEC error handling package.)
793 C
794 C Form of Call Function
795 C CALL XSETUN(LUN) set the logical unit number, LUN, for
796 C output of messages from DLSODA, if
797 C the default is not desired.
798 C The default value of LUN is 6.
799 C
800 C CALL XSETF(MFLAG) set a flag to control the printing of
801 C messages by DLSODA.
802 C MFLAG = 0 means do not print. (Danger:
803 C This risks losing valuable information.)
804 C MFLAG = 1 means print (the default).
805 C
806 C Either of the above calls may be made at
807 C any time and will take effect immediately.
808 C
809 C CALL DSRCMA(RSAV,ISAV,JOB) saves and restores the contents of
810 C the internal Common blocks used by
811 C DLSODA (see Part 3 below).
812 C RSAV must be a real array of length 240
813 C or more, and ISAV must be an integer
814 C array of length 46 or more.
815 C JOB=1 means save Common into RSAV/ISAV.
816 C JOB=2 means restore Common from RSAV/ISAV.
817 C DSRCMA is useful if one is
818 C interrupting a run and restarting
819 C later, or alternating between two or
820 C more problems solved with DLSODA.
821 C
822 C CALL DINTDY(,,,,,) provide derivatives of y, of various
823 C (see below) orders, at a specified point t, if
824 C desired. It may be called only after
825 C a successful return from DLSODA.
826 C
827 C The detailed instructions for using DINTDY are as follows.
828 C The form of the call is:
829 C
830 C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
831 C
832 C The input parameters are:
833 C
834 C T = value of independent variable where answers are desired
835 C (normally the same as the T last returned by DLSODA).
836 C For valid results, T must lie between TCUR - HU and TCUR.
837 C (See optional outputs for TCUR and HU.)
838 C K = integer order of the derivative desired. K must satisfy
839 C 0 .le. K .le. NQCUR, where NQCUR is the current order
840 C (see optional outputs). The capability corresponding
841 C to K = 0, i.e. computing y(T), is already provided
842 C by DLSODA directly. Since NQCUR .ge. 1, the first
843 C derivative dy/dt is always available with DINTDY.
844 C RWORK(21) = the base address of the history array YH.
845 C NYH = column length of YH, equal to the initial value of NEQ.
846 C
847 C The output parameters are:
848 C
849 C DKY = a real array of length NEQ containing the computed value
850 C of the K-th derivative of y(t).
851 C IFLAG = integer flag, returned as 0 if K and T were legal,
852 C -1 if K was illegal, and -2 if T was illegal.
853 C On an error return, a message is also written.
854 C-----------------------------------------------------------------------
855 C Part 3. Common Blocks.
856 C
857 C If DLSODA is to be used in an overlay situation, the user
858 C must declare, in the primary overlay, the variables in:
859 C (1) the call sequence to DLSODA, and
860 C (2) the two internal Common blocks
861 C /DLS001/ of length 255 (218 double precision words
862 C followed by 37 integer words),
863 C /DLSA01/ of length 31 (22 double precision words
864 C followed by 9 integer words).
865 C
866 C If DLSODA is used on a system in which the contents of internal
867 C Common blocks are not preserved between calls, the user should
868 C declare the above Common blocks in the calling program to insure
869 C that their contents are preserved.
870 C
871 C If the solution of a given problem by DLSODA is to be interrupted
872 C and then later continued, such as when restarting an interrupted run
873 C or alternating between two or more problems, the user should save,
874 C following the return from the last DLSODA call prior to the
875 C interruption, the contents of the call sequence variables and the
876 C internal Common blocks, and later restore these values before the
877 C next DLSODA call for that problem. To save and restore the Common
878 C blocks, use Subroutine DSRCMA (see Part 2 above).
879 C
880 C-----------------------------------------------------------------------
881 C Part 4. Optionally Replaceable Solver Routines.
882 C
883 C Below is a description of a routine in the DLSODA package which
884 C relates to the measurement of errors, and can be
885 C replaced by a user-supplied version, if desired. However, since such
886 C a replacement may have a major impact on performance, it should be
887 C done only when absolutely necessary, and only with great caution.
888 C (Note: The means by which the package version of a routine is
889 C superseded by the user's version may be system-dependent.)
890 C
891 C (a) DEWSET.
892 C The following subroutine is called just before each internal
893 C integration step, and sets the array of error weights, EWT, as
894 C described under ITOL/RTOL/ATOL above:
895 C Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
896 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODA call sequence,
897 C YCUR contains the current dependent variable vector, and
898 C EWT is the array of weights set by DEWSET.
899 C
900 C If the user supplies this subroutine, it must return in EWT(i)
901 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
902 C in y(i) to. The EWT array returned by DEWSET is passed to the
903 C DMNORM routine, and also used by DLSODA in the computation
904 C of the optional output IMXER, and the increments for difference
905 C quotient Jacobians.
906 C
907 C In the user-supplied version of DEWSET, it may be desirable to use
908 C the current values of derivatives of y. Derivatives up to order NQ
909 C are available from the history array YH, described above under
910 C optional outputs. In DEWSET, YH is identical to the YCUR array,
911 C extended to NQ + 1 columns with a column length of NYH and scale
912 C factors of H**j/factorial(j). On the first call for the problem,
913 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
914 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
915 C can be obtained by including in DEWSET the statements:
916 C DOUBLE PRECISION RLS
917 C COMMON /DLS001/ RLS(218),ILS(37)
918 C NQ = ILS(33)
919 C NST = ILS(34)
920 C H = RLS(212)
921 C Thus, for example, the current value of dy/dt can be obtained as
922 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
923 C unnecessary when NST = 0).
924 C-----------------------------------------------------------------------
925 C
926 C***REVISION HISTORY (YYYYMMDD)
927 C 19811102 DATE WRITTEN
928 C 19820126 Fixed bug in tests of work space lengths;
929 C minor corrections in main prologue and comments.
930 C 19870330 Major update: corrected comments throughout;
931 C removed TRET from Common; rewrote EWSET with 4 loops;
932 C fixed t test in INTDY; added Cray directives in STODA;
933 C in STODA, fixed DELP init. and logic around PJAC call;
934 C combined routines to save/restore Common;
935 C passed LEVEL = 0 in error message calls (except run abort).
936 C 19970225 Fixed lines setting JSTART = -2 in Subroutine LSODA.
937 C 20010425 Major update: convert source lines to upper case;
938 C added *DECK lines; changed from 1 to * in dummy dimensions;
939 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
940 C renamed routines for uniqueness across single/double prec.;
941 C converted intrinsic names to generic form;
942 C removed ILLIN and NTREP (data loaded) from Common;
943 C removed all 'own' variables from Common;
944 C changed error messages to quoted strings;
945 C replaced XERRWV/XERRWD with 1993 revised version;
946 C converted prologues, comments, error messages to mixed case;
947 C numerous corrections to prologues and internal comments.
948 C 20010507 Converted single precision source to double precision.
949 C 20010613 Revised excess accuracy test (to match rest of ODEPACK).
950 C 20010808 Fixed bug in DPRJA (matrix in DBNORM call).
951 C 20020502 Corrected declarations in descriptions of user routines.
952 C 20031105 Restored 'own' variables to Common blocks, to enable
953 C interrupt/restart feature.
954 C 20031112 Added SAVE statements for data-loaded constants.
955 C
956 C-----------------------------------------------------------------------
957 C Other routines in the DLSODA package.
958 C
959 C In addition to Subroutine DLSODA, the DLSODA package includes the
960 C following subroutines and function routines:
961 C DINTDY computes an interpolated value of the y vector at t = TOUT.
962 C DSTODA is the core integrator, which does one step of the
963 C integration and the associated error control.
964 C DCFODE sets all method coefficients and test constants.
965 C DPRJA computes and preprocesses the Jacobian matrix J = df/dy
966 C and the Newton iteration matrix P = I - h*l0*J.
967 C DSOLSY manages solution of linear system in chord iteration.
968 C DEWSET sets the error weight vector EWT before each step.
969 C DMNORM computes the weighted max-norm of a vector.
970 C DFNORM computes the norm of a full matrix consistent with the
971 C weighted max-norm on vectors.
972 C DBNORM computes the norm of a band matrix consistent with the
973 C weighted max-norm on vectors.
974 C DSRCMA is a user-callable routine to save and restore
975 C the contents of the internal Common blocks.
976 C DGEFA and DGESL are routines from LINPACK for solving full
977 C systems of linear algebraic equations.
978 C DGBFA and DGBSL are routines from LINPACK for solving banded
979 C linear systems.
980 C DUMACH computes the unit roundoff in a machine-independent manner.
981 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
982 C error messages and warnings. XERRWD is machine-dependent.
983 C Note: DMNORM, DFNORM, DBNORM, DUMACH, IXSAV, and IUMACH are
984 C function routines. All the others are subroutines.
985 C
986 C-----------------------------------------------------------------------
1003 LOGICAL IHIT
1004 CHARACTER*60 MSG
1006 C-----------------------------------------------------------------------
1007 C The following two internal Common blocks contain
1008 C (a) variables which are local to any subroutine but whose values must
1009 C be preserved between calls to the routine ("own" variables), and
1010 C (b) variables which are communicated between subroutines.
1011 C The block DLS001 is declared in subroutines DLSODA, DINTDY, DSTODA,
1012 C DPRJA, and DSOLSY.
1013 C The block DLSA01 is declared in subroutines DLSODA, DSTODA, and DPRJA.
1014 C Groups of variables are replaced by dummy arrays in the Common
1015 C declarations in routines where those variables are not used.
1016 C-----------------------------------------------------------------------
1023 C
1026 C
1028 C-----------------------------------------------------------------------
1029 C Block A.
1030 C This code block is executed on every call.
1031 C It tests ISTATE and ITASK for legality and branches appropriately.
1032 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1033 C not yet been done, an error return occurs.
1034 C If ISTATE = 1 and TOUT = T, return immediately.
1035 C-----------------------------------------------------------------------
1044 C-----------------------------------------------------------------------
1045 C Block B.
1046 C The next code block is executed for the initial call (ISTATE = 1),
1047 C or for a continuation call with parameter changes (ISTATE = 3).
1048 C It contains checking of all inputs and various initializations.
1049 C
1050 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1051 C JT, ML, and MU.
1052 C-----------------------------------------------------------------------
1060 JTYP = JT
1067 C Next process and check the optional inputs. --------------------------
1070 MXSTEP = MXSTP0
1071 MXHNIL = MXHNL0
1104 C-----------------------------------------------------------------------
1105 C Set work array pointers and check lengths LRW and LIW.
1106 C If ISTATE = 1, METH is initialized to 1 here to facilitate the
1107 C checking of work space lengths.
1108 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1109 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1110 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1111 C If the lengths provided are insufficient for the current method,
1112 C an error return occurs. This is treated as illegal input on the
1113 C first call, but as a problem interruption with ISTATE = -7 on a
1114 C continuation call. If the lengths are sufficient for the current
1115 C method but not for both methods, a warning message is sent.
1116 C-----------------------------------------------------------------------
1126 LEN1C = LEN1N
1165 C Check RTOL and ATOL for legality. ------------------------------------
1175 C If ISTATE = 3, set flag to signal parameter changes to DSTODA. -------
1178 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1185 C-----------------------------------------------------------------------
1186 C Block C.
1187 C The next block is for the initial call only (ISTATE = 1).
1188 C It contains all remaining initializations, the initial call to F,
1189 C and the calculation of the initial step size.
1190 C The error weights in EWT are inverted after being loaded.
1191 C-----------------------------------------------------------------------
1193 TN = T
1194 TSW = T
1195 MAXORD = MXORDN
1214 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1218 C Load the initial value vector in YH. ---------------------------------
1221 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1228 C-----------------------------------------------------------------------
1229 C The coding below computes the step size, H0, to be attempted on the
1230 C first step, unless the user has supplied a value for this.
1231 C First check that TOUT - T differs significantly from zero.
1232 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1233 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1234 C so as to be between 100*UROUND and 1.0E-3.
1235 C Then the computed value H0 is given by:
1236 C
1237 C H0**(-2) = 1./(TOL * w0**2) + TOL * (norm(F))**2
1238 C
1239 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1240 C F = the initial value of the vector f(t,y), and
1241 C norm() = the weighted vector norm used throughout, given by
1242 C the DMNORM function routine, and weighted by the
1243 C tolerances initially loaded into the EWT array.
1244 C The sign of H0 is inferred from the initial values of TOUT and T.
1245 C ABS(H0) is made .le. ABS(TOUT-T) in any case.
1246 C-----------------------------------------------------------------------
1269 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1272 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1273 H = H0
1277 C-----------------------------------------------------------------------
1278 C Block D.
1279 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1280 C and is to check stop conditions before taking a step.
1281 C-----------------------------------------------------------------------
1287 T = TOUT
1292 T = TN
1300 T = TOUT
1312 C-----------------------------------------------------------------------
1313 C Block E.
1314 C The next block is normally executed for all calls and contains
1315 C the call to the one-step core integrator DSTODA.
1316 C
1317 C This is a looping point for the integration steps.
1318 C
1319 C First check for too many steps being taken, update EWT (if not at
1320 C start of problem), check for too much accuracy being requested, and
1321 C check for H below the roundoff level in T.
1322 C-----------------------------------------------------------------------
1352 C-----------------------------------------------------------------------
1353 C CALL DSTODA(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPRJA,DSOLSY)
1354 C-----------------------------------------------------------------------
1360 C-----------------------------------------------------------------------
1361 C Block F.
1362 C The following block handles the case of a successful return from the
1363 C core integrator (KFLAG = 0).
1364 C If a method switch was just made, record TSW, reset MAXORD,
1365 C set JSTART to -1 to signal DSTODA to complete the switch,
1366 C and do extra printing of data if IXPR = 1.
1367 C Then, in any case, check for stop conditions.
1368 C-----------------------------------------------------------------------
1371 TSW = TN
1372 MAXORD = MXORDN
1382 ENDIF
1386 ENDIF
1390 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1393 T = TOUT
1395 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1398 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1401 T = TOUT
1411 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1414 C-----------------------------------------------------------------------
1415 C Block G.
1416 C The following block handles all successful returns from DLSODA.
1417 C If ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1418 C ISTATE is set to 2, and the optional outputs are loaded into the
1419 C work arrays before returning.
1420 C-----------------------------------------------------------------------
1423 T = TN
1438 RETURN
1439 C-----------------------------------------------------------------------
1440 C Block H.
1441 C The following block handles all unsuccessful returns other than
1442 C those for illegal input. First the error message routine is called.
1443 C If there was an error test or convergence test failure, IMXER is set.
1444 C Then Y is loaded from YH and T is set to TN.
1445 C The optional outputs are loaded into the work arrays before returning.
1446 C-----------------------------------------------------------------------
1447 C The maximum number of steps was taken before reaching TOUT. ----------
1454 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1460 C Too much accuracy requested for machine precision. -------------------
1468 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1475 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1484 C RWORK length too small to proceed. -----------------------------------
1491 C IWORK length too small to proceed. -----------------------------------
1498 C Compute IMXER if relevant. -------------------------------------------
1504 BIG = SIZE
1505 IMXER = I
1508 C Set Y vector, T, and optional outputs. -------------------------------
1511 T = TN
1523 RETURN
1524 C-----------------------------------------------------------------------
1525 C Block I.
1526 C The following block handles all error returns due to illegal input
1527 C (ISTATE = -3), as detected before calling the core integrator.
1528 C First the error message routine is called. If the illegal input
1529 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1530 C-----------------------------------------------------------------------
1624 C
1626 RETURN
1627 C
1630 RETURN
1631 C----------------------- End of Subroutine DLSODA ----------------------
1632 END