GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
slsode.f
Go to the documentation of this file.
1 *DECK SLSODE
2  SUBROUTINE slsode (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
3  1 istate, iopt, rwork, lrw, iwork, liw, jac, mf)
4  EXTERNAL f, jac
5  INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, MF
6  REAL Y, T, TOUT, RTOL, ATOL, RWORK
7  dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
8 C***BEGIN PROLOGUE SLSODE
9 C***PURPOSE Livermore Solver for Ordinary Differential Equations.
10 C SLSODE 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 SINGLE 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 SLSODE.
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 SLSODE 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 SLSODE package
35 C - internal COMMON block
36 C
37 C *Usage:
38 C Communication between the user and the SLSODE 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 REAL Y(NEQ), T, TOUT, RTOL, ATOL(ntol), RWORK(LRW)
61 C
62 C CALL SLSODE (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 REAL 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 SLSODE 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, as TOUT was equal to T.
130 C 2 SLSODE 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 REAL 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 SLSODE 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 SLSODE 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 REAL 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 SLSODE 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, SLSODE 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 SLSODE 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 REAL 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 SLSODE 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 SLSODE.
279 C
280 C Before the first call to SLSODE, 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 SLSODE 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 SLSODE. 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 SLSODE, 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 REAL ATOL(3), RTOL, RWORK(58), T, TOUT, Y(3)
308 C NEQ = 3
309 C Y(1) = 1.
310 C Y(2) = 0.
311 C Y(3) = 0.
312 C T = 0.
313 C TOUT = .4
314 C ITOL = 2
315 C RTOL = 1.E-4
316 C ATOL(1) = 1.E-6
317 C ATOL(2) = 1.E-10
318 C ATOL(3) = 1.E-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 SLSODE (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 =',E12.4,' y =',3E14.6)
330 C IF (ISTATE .LT. 0) GO TO 80
331 C 40 TOUT = TOUT*10.
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 REAL T, Y(3), YDOT(3)
343 C YDOT(1) = -.04*Y(1) + 1.E4*Y(2)*Y(3)
344 C YDOT(3) = 3.E7*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 REAL T, Y(3), PD(NRPD,3)
352 C PD(1,1) = -.04
353 C PD(1,2) = 1.E4*Y(3)
354 C PD(1,3) = 1.E4*Y(2)
355 C PD(2,1) = .04
356 C PD(2,3) = -PD(1,3)
357 C PD(3,2) = 6.E7*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 SLSODE 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 SLSODE, some compilers may object
392 C to a call to SLSODE 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 SLSODE consists of four parts:
407 C
408 C 1. The call sequence to subroutine SLSODE, 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 SLSODE 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 SLSODE 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 SLSODE 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 REAL 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 SLSODE, 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 SINTDY
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 SLSODE 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 SLSODE 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 (single 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 SLSODE
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 SLSODE 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 REAL 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 SLSODE.
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 SLSODE, the variables listed
880 C below are quantities related to the performance of SLSODE 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 SLSODE, 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 SLSODE.
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 SLSODE.
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 SLSODE, 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 SLSODE. 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 SSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the
970 C internal COMMON blocks used by SLSODE
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. SSRCOM 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 SLSODE.
980 C CALL SINTDY(,,,,,) 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 SLSODE. Detailed
984 C instructions follow.
985 C
986 C Detailed instructions for using SINTDY
987 C --------------------------------------
988 C The form of the CALL is:
989 C
990 C CALL SINTDY (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 SLSODE). 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 SLSODE directly. Since
1004 C NQCUR >= 1, the first derivative dy/dt is always
1005 C available with SINTDY.
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 SLSODE 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 SLSODE,
1024 C (2) the internal COMMON block /SLS001/, of length 255
1025 C (218 single precision words followed by 37 integer words).
1026 C
1027 C If SLSODE 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 SLSODE 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 SLSODE 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 SLSODE 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 SSRCOM (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 SLSODE 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 SEWSET
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 SEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1063 C
1064 C where NEQ, ITOL, RTOL, and ATOL are as in the SLSODE call
1065 C sequence, YCUR contains the current dependent variable vector,
1066 C and EWT is the array of weights set by SEWSET.
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 SEWSET is passed to the
1071 C SVNORM routine (see below), and also used by SLSODE 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 SEWSET, 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 SEWSET, 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 REAL RLS
1086 C COMMON /SLS001/ 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 SVNORM
1095 C ------
1096 C SVNORM is a real function routine which computes the weighted
1097 C root-mean-square norm of a vector v:
1098 C
1099 C d = SVNORM (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 SVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where
1108 C EWT is as set by subroutine SEWSET.
1109 C
1110 C If the user supplies this function, it should return a nonnegative
1111 C value of SVNORM suitable for use in the error control in SLSODE.
1112 C None of the arguments should be altered by SVNORM. For example, a
1113 C user-supplied SVNORM 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 SEWSET, SINTDY, R1MACH, SSTODE, SVNORM, XERRWD
1119 C***COMMON BLOCKS SLS001
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 XERRWV 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 XERRWV; 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 XERRWV 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 R1MACH to RUMACH. (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 XERRWV 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 SLSODE etc. (ACH)
1166 C 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH)
1167 C 20010412 Removed all 'own' variables from Common block /SLS001/
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 /SLS001/, to
1171 C enable interrupt/restart feature. (ACH)
1172 C 20031112 Added SAVE statements for data-loaded constants.
1173 C
1174 C*** END PROLOGUE SLSODE
1175 C
1176 C*Internal Notes:
1177 C
1178 C Other Routines in the SLSODE Package.
1179 C
1180 C In addition to Subroutine SLSODE, the SLSODE package includes the
1181 C following subroutines and function routines:
1182 C SINTDY computes an interpolated value of the y vector at t = TOUT.
1183 C SSTODE is the core integrator, which does one step of the
1184 C integration and the associated error control.
1185 C SCFODE sets all method coefficients and test constants.
1186 C SPREPJ computes and preprocesses the Jacobian matrix J = df/dy
1187 C and the Newton iteration matrix P = I - h*l0*J.
1188 C SSOLSY manages solution of linear system in chord iteration.
1189 C SEWSET sets the error weight vector EWT before each step.
1190 C SVNORM computes the weighted R.M.S. norm of a vector.
1191 C SSRCOM is a user-callable routine to save and restore
1192 C the contents of the internal Common block.
1193 C DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL
1194 C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS.
1195 C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED
1196 C LINEAR SYSTEMS.
1197 C R1MACH 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: SVNORM, R1MACH, IXSAV, and IUMACH are function routines.
1201 C All the others are subroutines.
1202 C
1203 C**End
1204 C
1205 C Declare externals.
1206  EXTERNAL sprepj, ssolsy
1207  REAL R1MACH, SVNORM
1208 C
1209 C Declare all other variables.
1210  INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH,
1211  1 ialth, ipup, lmax, meo, nqnyh, nslp,
1212  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1213  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1214  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1215  INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0,
1216  1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
1217  REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
1218  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1219  REAL ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI,
1220  1 tcrit, tdist, tnext, tol, tolsf, tp, SIZE, sum, w0
1221  dimension mord(2)
1222  LOGICAL IHIT
1223  CHARACTER*80 MSG
1224  SAVE mord, mxstp0, mxhnl0
1225 C-----------------------------------------------------------------------
1226 C The following internal Common block contains
1227 C (a) variables which are local to any subroutine but whose values must
1228 C be preserved between calls to the routine ("own" variables), and
1229 C (b) variables which are communicated between subroutines.
1230 C The block SLS001 is declared in subroutines SLSODE, SINTDY, SSTODE,
1231 C SPREPJ, and SSOLSY.
1232 C Groups of variables are replaced by dummy arrays in the Common
1233 C declarations in routines where those variables are not used.
1234 C-----------------------------------------------------------------------
1235  COMMON /sls001/ conit, crate, el(13), elco(13,12),
1236  1 hold, rmax, tesco(3,12),
1237  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1238  2 init, mxstep, mxhnil, nhnil, nslast, nyh,
1239  3 ialth, ipup, lmax, meo, nqnyh, nslp,
1240  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1241  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1242  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1243 C
1244  DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1245 C-----------------------------------------------------------------------
1246 C Block A.
1247 C This code block is executed on every call.
1248 C It tests ISTATE and ITASK for legality and branches appropriately.
1249 C If ISTATE .GT. 1 but the flag INIT shows that initialization has
1250 C not yet been done, an error return occurs.
1251 C If ISTATE = 1 and TOUT = T, return immediately.
1252 C-----------------------------------------------------------------------
1253 C
1254 C***FIRST EXECUTABLE STATEMENT SLSODE
1255  IF (istate .LT. 1 .OR. istate .GT. 3) go to 601
1256  IF (itask .LT. 1 .OR. itask .GT. 5) go to 602
1257  IF (istate .EQ. 1) go to 10
1258  IF (init .EQ. 0) go to 603
1259  IF (istate .EQ. 2) go to 200
1260  go to 20
1261  10 init = 0
1262  IF (tout .EQ. t) RETURN
1263 C-----------------------------------------------------------------------
1264 C Block B.
1265 C The next code block is executed for the initial call (ISTATE = 1),
1266 C or for a continuation call with parameter changes (ISTATE = 3).
1267 C It contains checking of all inputs and various initializations.
1268 C
1269 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1270 C MF, ML, and MU.
1271 C-----------------------------------------------------------------------
1272  20 IF (neq(1) .LE. 0) go to 604
1273  IF (istate .EQ. 1) go to 25
1274  IF (neq(1) .GT. n) go to 605
1275  25 n = neq(1)
1276  IF (itol .LT. 1 .OR. itol .GT. 4) go to 606
1277  IF (iopt .LT. 0 .OR. iopt .GT. 1) go to 607
1278  meth = mf/10
1279  miter = mf - 10*meth
1280  IF (meth .LT. 1 .OR. meth .GT. 2) go to 608
1281  IF (miter .LT. 0 .OR. miter .GT. 5) go to 608
1282  IF (miter .LE. 3) go to 30
1283  ml = iwork(1)
1284  mu = iwork(2)
1285  IF (ml .LT. 0 .OR. ml .GE. n) go to 609
1286  IF (mu .LT. 0 .OR. mu .GE. n) go to 610
1287  30 CONTINUE
1288 C Next process and check the optional inputs. --------------------------
1289  IF (iopt .EQ. 1) go to 40
1290  maxord = mord(meth)
1291  mxstep = mxstp0
1292  mxhnil = mxhnl0
1293  IF (istate .EQ. 1) h0 = 0.0e0
1294  hmxi = 0.0e0
1295  hmin = 0.0e0
1296  go to 60
1297  40 maxord = iwork(5)
1298  IF (maxord .LT. 0) go to 611
1299  IF (maxord .EQ. 0) maxord = 100
1300  maxord = min(maxord,mord(meth))
1301  mxstep = iwork(6)
1302  IF (mxstep .LT. 0) go to 612
1303  IF (mxstep .EQ. 0) mxstep = mxstp0
1304  mxhnil = iwork(7)
1305  IF (mxhnil .LT. 0) go to 613
1306  IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1307  IF (istate .NE. 1) go to 50
1308  h0 = rwork(5)
1309  IF ((tout - t)*h0 .LT. 0.0e0) go to 614
1310  50 hmax = rwork(6)
1311  IF (hmax .LT. 0.0e0) go to 615
1312  hmxi = 0.0e0
1313  IF (hmax .GT. 0.0e0) hmxi = 1.0e0/hmax
1314  hmin = rwork(7)
1315  IF (hmin .LT. 0.0e0) go to 616
1316 C-----------------------------------------------------------------------
1317 C Set work array pointers and check lengths LRW and LIW.
1318 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1319 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1320 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1321 C-----------------------------------------------------------------------
1322  60 lyh = 21
1323  IF (istate .EQ. 1) nyh = n
1324  lwm = lyh + (maxord + 1)*nyh
1325  IF (miter .EQ. 0) lenwm = 0
1326  IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1327  IF (miter .EQ. 3) lenwm = n + 2
1328  IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1329  lewt = lwm + lenwm
1330  lsavf = lewt + n
1331  lacor = lsavf + n
1332  lenrw = lacor + n - 1
1333  iwork(17) = lenrw
1334  liwm = 1
1335  leniw = 20 + n
1336  IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1337  iwork(18) = leniw
1338  IF (lenrw .GT. lrw) go to 617
1339  IF (leniw .GT. liw) go to 618
1340 C Check RTOL and ATOL for legality. ------------------------------------
1341  rtoli = rtol(1)
1342  atoli = atol(1)
1343  DO 70 i = 1,n
1344  IF (itol .GE. 3) rtoli = rtol(i)
1345  IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1346  IF (rtoli .LT. 0.0e0) go to 619
1347  IF (atoli .LT. 0.0e0) go to 620
1348  70 CONTINUE
1349  IF (istate .EQ. 1) go to 100
1350 C If ISTATE = 3, set flag to signal parameter changes to SSTODE. -------
1351  jstart = -1
1352  IF (nq .LE. maxord) go to 90
1353 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1354  DO 80 i = 1,n
1355  80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1356 C Reload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
1357  90 IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1358  IF (n .EQ. nyh) go to 200
1359 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1360  i1 = lyh + l*nyh
1361  i2 = lyh + (maxord + 1)*nyh - 1
1362  IF (i1 .GT. i2) go to 200
1363  DO 95 i = i1,i2
1364  95 rwork(i) = 0.0e0
1365  go to 200
1366 C-----------------------------------------------------------------------
1367 C Block C.
1368 C The next block is for the initial call only (ISTATE = 1).
1369 C It contains all remaining initializations, the initial call to F,
1370 C and the calculation of the initial step size.
1371 C The error weights in EWT are inverted after being loaded.
1372 C-----------------------------------------------------------------------
1373  100 uround = r1mach(4)
1374  tn = t
1375  IF (itask .NE. 4 .AND. itask .NE. 5) go to 110
1376  tcrit = rwork(1)
1377  IF ((tcrit - tout)*(tout - t) .LT. 0.0e0) go to 625
1378  IF (h0 .NE. 0.0e0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0e0)
1379  1 h0 = tcrit - t
1380  110 jstart = 0
1381  IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1382  nhnil = 0
1383  nst = 0
1384  nje = 0
1385  nslast = 0
1386  hu = 0.0e0
1387  nqu = 0
1388  ccmax = 0.3e0
1389  maxcor = 3
1390  msbp = 20
1391  mxncf = 10
1392 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1393  lf0 = lyh + nyh
1394  CALL f(neq, t, y, rwork(lf0))
1395  nfe = 1
1396 C Load the initial value vector in YH. ---------------------------------
1397  DO 115 i = 1,n
1398  115 rwork(i+lyh-1) = y(i)
1399 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1400  nq = 1
1401  h = 1.0e0
1402  CALL sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1403  DO 120 i = 1,n
1404  IF (rwork(i+lewt-1) .LE. 0.0e0) go to 621
1405  120 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1406 C-----------------------------------------------------------------------
1407 C The coding below computes the step size, H0, to be attempted on the
1408 C first step, unless the user has supplied a value for this.
1409 C First check that TOUT - T differs significantly from zero.
1410 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(I))
1411 C if this is positive, or MAX(ATOL(I)/ABS(Y(I))) otherwise, adjusted
1412 C so as to be between 100*UROUND and 1.0E-3.
1413 C Then the computed value H0 is given by..
1414 C NEQ
1415 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 )
1416 C 1
1417 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1418 C f(i) = i-th component of initial value of f,
1419 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1420 C The sign of H0 is inferred from the initial values of TOUT and T.
1421 C-----------------------------------------------------------------------
1422  IF (h0 .NE. 0.0e0) go to 180
1423  tdist = abs(tout - t)
1424  w0 = max(abs(t),abs(tout))
1425  IF (tdist .LT. 2.0e0*uround*w0) go to 622
1426  tol = rtol(1)
1427  IF (itol .LE. 2) go to 140
1428  DO 130 i = 1,n
1429  130 tol = max(tol,rtol(i))
1430  140 IF (tol .GT. 0.0e0) go to 160
1431  atoli = atol(1)
1432  DO 150 i = 1,n
1433  IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1434  ayi = abs(y(i))
1435  IF (ayi .NE. 0.0e0) tol = max(tol,atoli/ayi)
1436  150 CONTINUE
1437  160 tol = max(tol,100.0e0*uround)
1438  tol = min(tol,0.001e0)
1439  sum = svnorm(n, rwork(lf0), rwork(lewt))
1440  sum = 1.0e0/(tol*w0*w0) + tol*sum**2
1441  h0 = 1.0e0/sqrt(sum)
1442  h0 = min(h0,tdist)
1443  h0 = sign(h0,tout-t)
1444 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1445  180 rh = abs(h0)*hmxi
1446  IF (rh .GT. 1.0e0) h0 = h0/rh
1447 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1448  h = h0
1449  DO 190 i = 1,n
1450  190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1451  go to 270
1452 C-----------------------------------------------------------------------
1453 C Block D.
1454 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1455 C and is to check stop conditions before taking a step.
1456 C-----------------------------------------------------------------------
1457  200 nslast = nst
1458  go to(210, 250, 220, 230, 240), itask
1459  210 IF ((tn - tout)*h .LT. 0.0e0) go to 250
1460  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1461  IF (iflag .NE. 0) go to 627
1462  t = tout
1463  go to 420
1464  220 tp = tn - hu*(1.0e0 + 100.0e0*uround)
1465  IF ((tp - tout)*h .GT. 0.0e0) go to 623
1466  IF ((tn - tout)*h .LT. 0.0e0) go to 250
1467  go to 400
1468  230 tcrit = rwork(1)
1469  IF ((tn - tcrit)*h .GT. 0.0e0) go to 624
1470  IF ((tcrit - tout)*h .LT. 0.0e0) go to 625
1471  IF ((tn - tout)*h .LT. 0.0e0) go to 245
1472  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1473  IF (iflag .NE. 0) go to 627
1474  t = tout
1475  go to 420
1476  240 tcrit = rwork(1)
1477  IF ((tn - tcrit)*h .GT. 0.0e0) go to 624
1478  245 hmx = abs(tn) + abs(h)
1479  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1480  IF (ihit) go to 400
1481  tnext = tn + h*(1.0e0 + 4.0e0*uround)
1482  IF ((tnext - tcrit)*h .LE. 0.0e0) go to 250
1483  h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1484  IF (istate .EQ. 2) jstart = -2
1485 C-----------------------------------------------------------------------
1486 C Block E.
1487 C The next block is normally executed for all calls and contains
1488 C the call to the one-step core integrator SSTODE.
1489 C
1490 C This is a looping point for the integration steps.
1491 C
1492 C First check for too many steps being taken, update EWT (if not at
1493 C start of problem), check for too much accuracy being requested, and
1494 C check for H below the roundoff level in T.
1495 C-----------------------------------------------------------------------
1496  250 CONTINUE
1497  IF ((nst-nslast) .GE. mxstep) go to 500
1498  CALL sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1499  DO 260 i = 1,n
1500  IF (rwork(i+lewt-1) .LE. 0.0e0) go to 510
1501  260 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1502  270 tolsf = uround*svnorm(n, rwork(lyh), rwork(lewt))
1503  IF (tolsf .LE. 1.0e0) go to 280
1504  tolsf = tolsf*2.0e0
1505  IF (nst .EQ. 0) go to 626
1506  go to 520
1507  280 IF ((tn + h) .NE. tn) go to 290
1508  nhnil = nhnil + 1
1509  IF (nhnil .GT. mxhnil) go to 290
1510  CALL xerrwd('SLSODE- Warning..internal T (=R1) and H (=R2) are',
1511  1 50, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1512  CALL xerrwd(
1513  1 ' such that in the machine, T + H = T on the next step ',
1514  1 60, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1515  CALL xerrwd(' (H = step size). Solver will continue anyway',
1516  1 50, 101, 0, 0, 0, 0, 2, tn, h)
1517  IF (nhnil .LT. mxhnil) go to 290
1518  CALL xerrwd('SLSODE- Above warning has been issued I1 times. ',
1519  1 50, 102, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1520  CALL xerrwd(' It will not be issued again for this problem',
1521  1 50, 102, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1522  290 CONTINUE
1523 C-----------------------------------------------------------------------
1524 C CALL SSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,SPREPJ,SSOLSY)
1525 C-----------------------------------------------------------------------
1526  CALL sstode(neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1527  1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1528  2 f, jac, sprepj, ssolsy)
1529  kgo = 1 - kflag
1530  go to(300, 530, 540), kgo
1531 C-----------------------------------------------------------------------
1532 C Block F.
1533 C The following block handles the case of a successful return from the
1534 C core integrator (KFLAG = 0). Test for stop conditions.
1535 C-----------------------------------------------------------------------
1536  300 init = 1
1537  go to(310, 400, 330, 340, 350), itask
1538 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1539  310 IF ((tn - tout)*h .LT. 0.0e0) go to 250
1540  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1541  t = tout
1542  go to 420
1543 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1544  330 IF ((tn - tout)*h .GE. 0.0e0) go to 400
1545  go to 250
1546 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1547  340 IF ((tn - tout)*h .LT. 0.0e0) go to 345
1548  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1549  t = tout
1550  go to 420
1551  345 hmx = abs(tn) + abs(h)
1552  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1553  IF (ihit) go to 400
1554  tnext = tn + h*(1.0e0 + 4.0e0*uround)
1555  IF ((tnext - tcrit)*h .LE. 0.0e0) go to 250
1556  h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1557  jstart = -2
1558  go to 250
1559 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1560  350 hmx = abs(tn) + abs(h)
1561  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1562 C-----------------------------------------------------------------------
1563 C Block G.
1564 C The following block handles all successful returns from SLSODE.
1565 C If ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
1566 C ISTATE is set to 2, and the optional outputs are loaded into the
1567 C work arrays before returning.
1568 C-----------------------------------------------------------------------
1569  400 DO 410 i = 1,n
1570  410 y(i) = rwork(i+lyh-1)
1571  t = tn
1572  IF (itask .NE. 4 .AND. itask .NE. 5) go to 420
1573  IF (ihit) t = tcrit
1574  420 istate = 2
1575  rwork(11) = hu
1576  rwork(12) = h
1577  rwork(13) = tn
1578  iwork(11) = nst
1579  iwork(12) = nfe
1580  iwork(13) = nje
1581  iwork(14) = nqu
1582  iwork(15) = nq
1583  RETURN
1584 C-----------------------------------------------------------------------
1585 C Block H.
1586 C The following block handles all unsuccessful returns other than
1587 C those for illegal input. First the error message routine is called.
1588 C If there was an error test or convergence test failure, IMXER is set.
1589 C Then Y is loaded from YH and T is set to TN. The optional outputs
1590 C are loaded into the work arrays before returning.
1591 C-----------------------------------------------------------------------
1592 C The maximum number of steps was taken before reaching TOUT. ----------
1593  500 CALL xerrwd('SLSODE- At current T (=R1), MXSTEP (=I1) steps ',
1594  1 50, 201, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1595  CALL xerrwd(' taken on this call before reaching TOUT ',
1596  1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0e0)
1597  istate = -1
1598  go to 580
1599 C EWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
1600  510 ewti = rwork(lewt+i-1)
1601  CALL xerrwd(.LE.'SLSODE- At T (=R1), EWT(I1) has become R2 0.',
1602  1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1603  istate = -6
1604  go to 580
1605 C Too much accuracy requested for machine precision. -------------------
1606  520 CALL xerrwd('SLSODE- At T (=R1), too much accuracy requested ',
1607  1 50, 203, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1608  CALL xerrwd(' for precision of machine.. see TOLSF (=R2) ',
1609  1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1610  rwork(14) = tolsf
1611  istate = -2
1612  go to 580
1613 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1614  530 CALL xerrwd('SLSODE- At T(=R1) and step size H(=R2), the error',
1615  1 50, 204, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1616  CALL xerrwd(' test failed repeatedly or with ABS(H) = HMIN',
1617  1 50, 204, 0, 0, 0, 0, 2, tn, h)
1618  istate = -4
1619  go to 560
1620 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1621  540 CALL xerrwd('SLSODE- At T (=R1) and step size H (=R2), the ',
1622  1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1623  CALL xerrwd(' corrector convergence failed repeatedly ',
1624  1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1625  CALL xerrwd(' or with ABS(H) = HMIN ',
1626  1 30, 205, 0, 0, 0, 0, 2, tn, h)
1627  istate = -5
1628 C Compute IMXER if relevant. -------------------------------------------
1629  560 big = 0.0e0
1630  imxer = 1
1631  DO 570 i = 1,n
1632  SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1633  IF (big .GE. size) go to 570
1634  big = SIZE
1635  imxer = i
1636  570 CONTINUE
1637  iwork(16) = imxer
1638 C Set Y vector, T, and optional outputs. -------------------------------
1639  580 DO 590 i = 1,n
1640  590 y(i) = rwork(i+lyh-1)
1641  t = tn
1642  rwork(11) = hu
1643  rwork(12) = h
1644  rwork(13) = tn
1645  iwork(11) = nst
1646  iwork(12) = nfe
1647  iwork(13) = nje
1648  iwork(14) = nqu
1649  iwork(15) = nq
1650  RETURN
1651 C-----------------------------------------------------------------------
1652 C Block I.
1653 C The following block handles all error returns due to illegal input
1654 C (ISTATE = -3), as detected before calling the core integrator.
1655 C First the error message routine is called. If the illegal input
1656 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1657 C-----------------------------------------------------------------------
1658  601 CALL xerrwd('SLSODE- ISTATE (=I1) illegal ',
1659  1 30, 1, 0, 1, istate, 0, 0, 0.0e0, 0.0e0)
1660  IF (istate .LT. 0) go to 800
1661  go to 700
1662  602 CALL xerrwd('SLSODE- ITASK (=I1) illegal ',
1663  1 30, 2, 0, 1, itask, 0, 0, 0.0e0, 0.0e0)
1664  go to 700
1665  603 CALL xerrwd(.GT.'SLSODE- ISTATE 1 but SLSODE not initialized ',
1666  1 50, 3, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1667  go to 700
1668  604 CALL xerrwd(.LT.'SLSODE- NEQ (=I1) 1 ',
1669  1 30, 4, 0, 1, neq(1), 0, 0, 0.0e0, 0.0e0)
1670  go to 700
1671  605 CALL xerrwd('SLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ',
1672  1 50, 5, 0, 2, n, neq(1), 0, 0.0e0, 0.0e0)
1673  go to 700
1674  606 CALL xerrwd('SLSODE- ITOL (=I1) illegal ',
1675  1 30, 6, 0, 1, itol, 0, 0, 0.0e0, 0.0e0)
1676  go to 700
1677  607 CALL xerrwd('SLSODE- IOPT (=I1) illegal ',
1678  1 30, 7, 0, 1, iopt, 0, 0, 0.0e0, 0.0e0)
1679  go to 700
1680  608 CALL xerrwd('SLSODE- MF (=I1) illegal ',
1681  1 30, 8, 0, 1, mf, 0, 0, 0.0e0, 0.0e0)
1682  go to 700
1683  609 CALL xerrwd(.LT..GE.'SLSODE- ML (=I1) illegal.. 0 or NEQ (=I2)',
1684  1 50, 9, 0, 2, ml, neq(1), 0, 0.0e0, 0.0e0)
1685  go to 700
1686  610 CALL xerrwd(.LT..GE.'SLSODE- MU (=I1) illegal.. 0 or NEQ (=I2)',
1687  1 50, 10, 0, 2, mu, neq(1), 0, 0.0e0, 0.0e0)
1688  go to 700
1689  611 CALL xerrwd(.LT.'SLSODE- MAXORD (=I1) 0 ',
1690  1 30, 11, 0, 1, maxord, 0, 0, 0.0e0, 0.0e0)
1691  go to 700
1692  612 CALL xerrwd(.LT.'SLSODE- MXSTEP (=I1) 0 ',
1693  1 30, 12, 0, 1, mxstep, 0, 0, 0.0e0, 0.0e0)
1694  go to 700
1695  613 CALL xerrwd(.LT.'SLSODE- MXHNIL (=I1) 0 ',
1696  1 30, 13, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1697  go to 700
1698  614 CALL xerrwd('SLSODE- TOUT (=R1) behind T (=R2) ',
1699  1 40, 14, 0, 0, 0, 0, 2, tout, t)
1700  CALL xerrwd(' Integration direction is given by H0 (=R1) ',
1701  1 50, 14, 0, 0, 0, 0, 1, h0, 0.0e0)
1702  go to 700
1703  615 CALL xerrwd(.LT.'SLSODE- HMAX (=R1) 0.0 ',
1704  1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0e0)
1705  go to 700
1706  616 CALL xerrwd(.LT.'SLSODE- HMIN (=R1) 0.0 ',
1707  1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0e0)
1708  go to 700
1709  617 CALL xerrwd(
1710  1 'SLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)',
1711  1 60, 17, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1712  go to 700
1713  618 CALL xerrwd(
1714  1 'SLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)',
1715  1 60, 18, 0, 2, leniw, liw, 0, 0.0e0, 0.0e0)
1716  go to 700
1717  619 CALL xerrwd(.LT.'SLSODE- RTOL(I1) is R1 0.0 ',
1718  1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0e0)
1719  go to 700
1720  620 CALL xerrwd(.LT.'SLSODE- ATOL(I1) is R1 0.0 ',
1721  1 40, 20, 0, 1, i, 0, 1, atoli, 0.0e0)
1722  go to 700
1723  621 ewti = rwork(lewt+i-1)
1724  CALL xerrwd(.LE.'SLSODE- EWT(I1) is R1 0.0 ',
1725  1 40, 21, 0, 1, i, 0, 1, ewti, 0.0e0)
1726  go to 700
1727  622 CALL xerrwd(
1728  1 'SLSODE- TOUT (=R1) too close to T(=R2) to start integration',
1729  1 60, 22, 0, 0, 0, 0, 2, tout, t)
1730  go to 700
1731  623 CALL xerrwd(
1732  1 'SLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ',
1733  1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1734  go to 700
1735  624 CALL xerrwd(
1736  1 'SLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ',
1737  1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1738  go to 700
1739  625 CALL xerrwd(
1740  1 'SLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ',
1741  1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1742  go to 700
1743  626 CALL xerrwd('SLSODE- At start of problem, too much accuracy ',
1744  1 50, 26, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1745  CALL xerrwd(
1746  1 ' requested for precision of machine.. See TOLSF (=R1) ',
1747  1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0e0)
1748  rwork(14) = tolsf
1749  go to 700
1750  627 CALL xerrwd('SLSODE- Trouble in SINTDY. ITASK = I1, TOUT = R1',
1751  1 50, 27, 0, 1, itask, 0, 1, tout, 0.0e0)
1752 C
1753  700 istate = -3
1754  RETURN
1755 C
1756  800 CALL xerrwd('SLSODE- Run aborted.. apparent infinite loop ',
1757  1 50, 303, 2, 0, 0, 0, 0, 0.0e0, 0.0e0)
1758  RETURN
1759 C----------------------- END OF SUBROUTINE SLSODE ----------------------
1760  END
subroutine ssolsy(WM, IWM, X, TEM)
Definition: ssolsy.f:1
subroutine sewset(N, ITOL, RTOL, ATOL, YCUR, EWT)
Definition: sewset.f:1
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:3
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
subroutine slsode(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
Definition: slsode.f:2
subroutine sintdy(T, K, YH, NYH, DKY, IFLAG)
Definition: sintdy.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
OCTAVE_EXPORT octave_value_list any number nd example oindent prints the prompt xample Pick a any number!nd example oindent and waits for the user to enter a value The string entered by the user is evaluated as an so it may be a literal a variable or any other valid Octave code The number of return their size
Definition: input.cc:871
subroutine sprepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
Definition: sprepj.f:1
subroutine sstode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
Definition: sstode.f:1
octave_value sqrt(void) const
Definition: ov.h:1388
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205