GNU Octave  4.0.0
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, D1MACH, 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 D1MACH 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, D1MACH, 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 D1MACH, SVNORM
1208 C
1209 C Declare all other variables.
1210  INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS,
1211  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1212  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1213  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1214  INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0,
1215  1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
1216  REAL ROWNS,
1217  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1218  REAL ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI,
1219  1 tcrit, tdist, tnext, tol, tolsf, tp, SIZE, sum, w0
1220  dimension mord(2)
1221  LOGICAL IHIT
1222  CHARACTER*80 MSG
1223  SAVE mord, mxstp0, mxhnl0
1224 C-----------------------------------------------------------------------
1225 C The following internal Common block contains
1226 C (a) variables which are local to any subroutine but whose values must
1227 C be preserved between calls to the routine ("own" variables), and
1228 C (b) variables which are communicated between subroutines.
1229 C The block SLS001 is declared in subroutines SLSODE, SINTDY, SSTODE,
1230 C SPREPJ, and SSOLSY.
1231 C Groups of variables are replaced by dummy arrays in the Common
1232 C declarations in routines where those variables are not used.
1233 C-----------------------------------------------------------------------
1234  COMMON /sls001/ rowns(209),
1235  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1236  2 init, mxstep, mxhnil, nhnil, nslast, nyh, iowns(6),
1237  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1238  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1239  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1240 C
1241  DATA mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1242 C-----------------------------------------------------------------------
1243 C Block A.
1244 C This code block is executed on every call.
1245 C It tests ISTATE and ITASK for legality and branches appropriately.
1246 C If ISTATE .GT. 1 but the flag INIT shows that initialization has
1247 C not yet been done, an error return occurs.
1248 C If ISTATE = 1 and TOUT = T, return immediately.
1249 C-----------------------------------------------------------------------
1250 C
1251 C***FIRST EXECUTABLE STATEMENT SLSODE
1252  IF (istate .LT. 1 .OR. istate .GT. 3) go to 601
1253  IF (itask .LT. 1 .OR. itask .GT. 5) go to 602
1254  IF (istate .EQ. 1) go to 10
1255  IF (init .EQ. 0) go to 603
1256  IF (istate .EQ. 2) go to 200
1257  go to 20
1258  10 init = 0
1259  IF (tout .EQ. t) RETURN
1260 C-----------------------------------------------------------------------
1261 C Block B.
1262 C The next code block is executed for the initial call (ISTATE = 1),
1263 C or for a continuation call with parameter changes (ISTATE = 3).
1264 C It contains checking of all inputs and various initializations.
1265 C
1266 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1267 C MF, ML, and MU.
1268 C-----------------------------------------------------------------------
1269  20 IF (neq(1) .LE. 0) go to 604
1270  IF (istate .EQ. 1) go to 25
1271  IF (neq(1) .GT. n) go to 605
1272  25 n = neq(1)
1273  IF (itol .LT. 1 .OR. itol .GT. 4) go to 606
1274  IF (iopt .LT. 0 .OR. iopt .GT. 1) go to 607
1275  meth = mf/10
1276  miter = mf - 10*meth
1277  IF (meth .LT. 1 .OR. meth .GT. 2) go to 608
1278  IF (miter .LT. 0 .OR. miter .GT. 5) go to 608
1279  IF (miter .LE. 3) go to 30
1280  ml = iwork(1)
1281  mu = iwork(2)
1282  IF (ml .LT. 0 .OR. ml .GE. n) go to 609
1283  IF (mu .LT. 0 .OR. mu .GE. n) go to 610
1284  30 CONTINUE
1285 C Next process and check the optional inputs. --------------------------
1286  IF (iopt .EQ. 1) go to 40
1287  maxord = mord(meth)
1288  mxstep = mxstp0
1289  mxhnil = mxhnl0
1290  IF (istate .EQ. 1) h0 = 0.0e0
1291  hmxi = 0.0e0
1292  hmin = 0.0e0
1293  go to 60
1294  40 maxord = iwork(5)
1295  IF (maxord .LT. 0) go to 611
1296  IF (maxord .EQ. 0) maxord = 100
1297  maxord = min(maxord,mord(meth))
1298  mxstep = iwork(6)
1299  IF (mxstep .LT. 0) go to 612
1300  IF (mxstep .EQ. 0) mxstep = mxstp0
1301  mxhnil = iwork(7)
1302  IF (mxhnil .LT. 0) go to 613
1303  IF (mxhnil .EQ. 0) mxhnil = mxhnl0
1304  IF (istate .NE. 1) go to 50
1305  h0 = rwork(5)
1306  IF ((tout - t)*h0 .LT. 0.0e0) go to 614
1307  50 hmax = rwork(6)
1308  IF (hmax .LT. 0.0e0) go to 615
1309  hmxi = 0.0e0
1310  IF (hmax .GT. 0.0e0) hmxi = 1.0e0/hmax
1311  hmin = rwork(7)
1312  IF (hmin .LT. 0.0e0) go to 616
1313 C-----------------------------------------------------------------------
1314 C Set work array pointers and check lengths LRW and LIW.
1315 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1316 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1317 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1318 C-----------------------------------------------------------------------
1319  60 lyh = 21
1320  IF (istate .EQ. 1) nyh = n
1321  lwm = lyh + (maxord + 1)*nyh
1322  IF (miter .EQ. 0) lenwm = 0
1323  IF (miter .EQ. 1 .OR. miter .EQ. 2) lenwm = n*n + 2
1324  IF (miter .EQ. 3) lenwm = n + 2
1325  IF (miter .GE. 4) lenwm = (2*ml + mu + 1)*n + 2
1326  lewt = lwm + lenwm
1327  lsavf = lewt + n
1328  lacor = lsavf + n
1329  lenrw = lacor + n - 1
1330  iwork(17) = lenrw
1331  liwm = 1
1332  leniw = 20 + n
1333  IF (miter .EQ. 0 .OR. miter .EQ. 3) leniw = 20
1334  iwork(18) = leniw
1335  IF (lenrw .GT. lrw) go to 617
1336  IF (leniw .GT. liw) go to 618
1337 C Check RTOL and ATOL for legality. ------------------------------------
1338  rtoli = rtol(1)
1339  atoli = atol(1)
1340  DO 70 i = 1,n
1341  IF (itol .GE. 3) rtoli = rtol(i)
1342  IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1343  IF (rtoli .LT. 0.0e0) go to 619
1344  IF (atoli .LT. 0.0e0) go to 620
1345  70 CONTINUE
1346  IF (istate .EQ. 1) go to 100
1347 C If ISTATE = 3, set flag to signal parameter changes to SSTODE. -------
1348  jstart = -1
1349  IF (nq .LE. maxord) go to 90
1350 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1351  DO 80 i = 1,n
1352  80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1353 C Reload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
1354  90 IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1355  IF (n .EQ. nyh) go to 200
1356 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1357  i1 = lyh + l*nyh
1358  i2 = lyh + (maxord + 1)*nyh - 1
1359  IF (i1 .GT. i2) go to 200
1360  DO 95 i = i1,i2
1361  95 rwork(i) = 0.0e0
1362  go to 200
1363 C-----------------------------------------------------------------------
1364 C Block C.
1365 C The next block is for the initial call only (ISTATE = 1).
1366 C It contains all remaining initializations, the initial call to F,
1367 C and the calculation of the initial step size.
1368 C The error weights in EWT are inverted after being loaded.
1369 C-----------------------------------------------------------------------
1370  100 uround = d1mach(4)
1371  tn = t
1372  IF (itask .NE. 4 .AND. itask .NE. 5) go to 110
1373  tcrit = rwork(1)
1374  IF ((tcrit - tout)*(tout - t) .LT. 0.0e0) go to 625
1375  IF (h0 .NE. 0.0e0 .AND. (t + h0 - tcrit)*h0 .GT. 0.0e0)
1376  1 h0 = tcrit - t
1377  110 jstart = 0
1378  IF (miter .GT. 0) rwork(lwm) = sqrt(uround)
1379  nhnil = 0
1380  nst = 0
1381  nje = 0
1382  nslast = 0
1383  hu = 0.0e0
1384  nqu = 0
1385  ccmax = 0.3e0
1386  maxcor = 3
1387  msbp = 20
1388  mxncf = 10
1389 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1390  lf0 = lyh + nyh
1391  CALL f(neq, t, y, rwork(lf0))
1392  nfe = 1
1393 C Load the initial value vector in YH. ---------------------------------
1394  DO 115 i = 1,n
1395  115 rwork(i+lyh-1) = y(i)
1396 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1397  nq = 1
1398  h = 1.0e0
1399  CALL sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1400  DO 120 i = 1,n
1401  IF (rwork(i+lewt-1) .LE. 0.0e0) go to 621
1402  120 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1403 C-----------------------------------------------------------------------
1404 C The coding below computes the step size, H0, to be attempted on the
1405 C first step, unless the user has supplied a value for this.
1406 C First check that TOUT - T differs significantly from zero.
1407 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(I))
1408 C if this is positive, or MAX(ATOL(I)/ABS(Y(I))) otherwise, adjusted
1409 C so as to be between 100*UROUND and 1.0E-3.
1410 C Then the computed value H0 is given by..
1411 C NEQ
1412 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 )
1413 C 1
1414 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1415 C f(i) = i-th component of initial value of f,
1416 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1417 C The sign of H0 is inferred from the initial values of TOUT and T.
1418 C-----------------------------------------------------------------------
1419  IF (h0 .NE. 0.0e0) go to 180
1420  tdist = abs(tout - t)
1421  w0 = max(abs(t),abs(tout))
1422  IF (tdist .LT. 2.0e0*uround*w0) go to 622
1423  tol = rtol(1)
1424  IF (itol .LE. 2) go to 140
1425  DO 130 i = 1,n
1426  130 tol = max(tol,rtol(i))
1427  140 IF (tol .GT. 0.0e0) go to 160
1428  atoli = atol(1)
1429  DO 150 i = 1,n
1430  IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
1431  ayi = abs(y(i))
1432  IF (ayi .NE. 0.0e0) tol = max(tol,atoli/ayi)
1433  150 CONTINUE
1434  160 tol = max(tol,100.0e0*uround)
1435  tol = min(tol,0.001e0)
1436  sum = svnorm(n, rwork(lf0), rwork(lewt))
1437  sum = 1.0e0/(tol*w0*w0) + tol*sum**2
1438  h0 = 1.0e0/sqrt(sum)
1439  h0 = min(h0,tdist)
1440  h0 = sign(h0,tout-t)
1441 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1442  180 rh = abs(h0)*hmxi
1443  IF (rh .GT. 1.0e0) h0 = h0/rh
1444 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1445  h = h0
1446  DO 190 i = 1,n
1447  190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1448  go to 270
1449 C-----------------------------------------------------------------------
1450 C Block D.
1451 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1452 C and is to check stop conditions before taking a step.
1453 C-----------------------------------------------------------------------
1454  200 nslast = nst
1455  go to(210, 250, 220, 230, 240), itask
1456  210 IF ((tn - tout)*h .LT. 0.0e0) go to 250
1457  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1458  IF (iflag .NE. 0) go to 627
1459  t = tout
1460  go to 420
1461  220 tp = tn - hu*(1.0e0 + 100.0e0*uround)
1462  IF ((tp - tout)*h .GT. 0.0e0) go to 623
1463  IF ((tn - tout)*h .LT. 0.0e0) go to 250
1464  go to 400
1465  230 tcrit = rwork(1)
1466  IF ((tn - tcrit)*h .GT. 0.0e0) go to 624
1467  IF ((tcrit - tout)*h .LT. 0.0e0) go to 625
1468  IF ((tn - tout)*h .LT. 0.0e0) go to 245
1469  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1470  IF (iflag .NE. 0) go to 627
1471  t = tout
1472  go to 420
1473  240 tcrit = rwork(1)
1474  IF ((tn - tcrit)*h .GT. 0.0e0) go to 624
1475  245 hmx = abs(tn) + abs(h)
1476  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1477  IF (ihit) go to 400
1478  tnext = tn + h*(1.0e0 + 4.0e0*uround)
1479  IF ((tnext - tcrit)*h .LE. 0.0e0) go to 250
1480  h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1481  IF (istate .EQ. 2) jstart = -2
1482 C-----------------------------------------------------------------------
1483 C Block E.
1484 C The next block is normally executed for all calls and contains
1485 C the call to the one-step core integrator SSTODE.
1486 C
1487 C This is a looping point for the integration steps.
1488 C
1489 C First check for too many steps being taken, update EWT (if not at
1490 C start of problem), check for too much accuracy being requested, and
1491 C check for H below the roundoff level in T.
1492 C-----------------------------------------------------------------------
1493  250 CONTINUE
1494  IF ((nst-nslast) .GE. mxstep) go to 500
1495  CALL sewset(n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1496  DO 260 i = 1,n
1497  IF (rwork(i+lewt-1) .LE. 0.0e0) go to 510
1498  260 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1499  270 tolsf = uround*svnorm(n, rwork(lyh), rwork(lewt))
1500  IF (tolsf .LE. 1.0e0) go to 280
1501  tolsf = tolsf*2.0e0
1502  IF (nst .EQ. 0) go to 626
1503  go to 520
1504  280 IF ((tn + h) .NE. tn) go to 290
1505  nhnil = nhnil + 1
1506  IF (nhnil .GT. mxhnil) go to 290
1507  CALL xerrwd('SLSODE- Warning..internal T (=R1) and H (=R2) are',
1508  1 50, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1509  CALL xerrwd(
1510  1 ' such that in the machine, T + H = T on the next step ',
1511  1 60, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1512  CALL xerrwd(' (H = step size). Solver will continue anyway',
1513  1 50, 101, 0, 0, 0, 0, 2, tn, h)
1514  IF (nhnil .LT. mxhnil) go to 290
1515  CALL xerrwd('SLSODE- Above warning has been issued I1 times. ',
1516  1 50, 102, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1517  CALL xerrwd(' It will not be issued again for this problem',
1518  1 50, 102, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1519  290 CONTINUE
1520 C-----------------------------------------------------------------------
1521 C CALL SSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,SPREPJ,SSOLSY)
1522 C-----------------------------------------------------------------------
1523  CALL sstode(neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1524  1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1525  2 f, jac, sprepj, ssolsy)
1526  kgo = 1 - kflag
1527  go to(300, 530, 540), kgo
1528 C-----------------------------------------------------------------------
1529 C Block F.
1530 C The following block handles the case of a successful return from the
1531 C core integrator (KFLAG = 0). Test for stop conditions.
1532 C-----------------------------------------------------------------------
1533  300 init = 1
1534  go to(310, 400, 330, 340, 350), itask
1535 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1536  310 IF ((tn - tout)*h .LT. 0.0e0) go to 250
1537  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1538  t = tout
1539  go to 420
1540 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1541  330 IF ((tn - tout)*h .GE. 0.0e0) go to 400
1542  go to 250
1543 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1544  340 IF ((tn - tout)*h .LT. 0.0e0) go to 345
1545  CALL sintdy(tout, 0, rwork(lyh), nyh, y, iflag)
1546  t = tout
1547  go to 420
1548  345 hmx = abs(tn) + abs(h)
1549  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1550  IF (ihit) go to 400
1551  tnext = tn + h*(1.0e0 + 4.0e0*uround)
1552  IF ((tnext - tcrit)*h .LE. 0.0e0) go to 250
1553  h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1554  jstart = -2
1555  go to 250
1556 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1557  350 hmx = abs(tn) + abs(h)
1558  ihit = abs(tn - tcrit) .LE. 100.0e0*uround*hmx
1559 C-----------------------------------------------------------------------
1560 C Block G.
1561 C The following block handles all successful returns from SLSODE.
1562 C If ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
1563 C ISTATE is set to 2, and the optional outputs are loaded into the
1564 C work arrays before returning.
1565 C-----------------------------------------------------------------------
1566  400 DO 410 i = 1,n
1567  410 y(i) = rwork(i+lyh-1)
1568  t = tn
1569  IF (itask .NE. 4 .AND. itask .NE. 5) go to 420
1570  IF (ihit) t = tcrit
1571  420 istate = 2
1572  rwork(11) = hu
1573  rwork(12) = h
1574  rwork(13) = tn
1575  iwork(11) = nst
1576  iwork(12) = nfe
1577  iwork(13) = nje
1578  iwork(14) = nqu
1579  iwork(15) = nq
1580  RETURN
1581 C-----------------------------------------------------------------------
1582 C Block H.
1583 C The following block handles all unsuccessful returns other than
1584 C those for illegal input. First the error message routine is called.
1585 C If there was an error test or convergence test failure, IMXER is set.
1586 C Then Y is loaded from YH and T is set to TN. The optional outputs
1587 C are loaded into the work arrays before returning.
1588 C-----------------------------------------------------------------------
1589 C The maximum number of steps was taken before reaching TOUT. ----------
1590  500 CALL xerrwd('SLSODE- At current T (=R1), MXSTEP (=I1) steps ',
1591  1 50, 201, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1592  CALL xerrwd(' taken on this call before reaching TOUT ',
1593  1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0e0)
1594  istate = -1
1595  go to 580
1596 C EWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
1597  510 ewti = rwork(lewt+i-1)
1598  CALL xerrwd(.LE.'SLSODE- At T (=R1), EWT(I1) has become R2 0.',
1599  1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1600  istate = -6
1601  go to 580
1602 C Too much accuracy requested for machine precision. -------------------
1603  520 CALL xerrwd('SLSODE- At T (=R1), too much accuracy requested ',
1604  1 50, 203, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1605  CALL xerrwd(' for precision of machine.. see TOLSF (=R2) ',
1606  1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1607  rwork(14) = tolsf
1608  istate = -2
1609  go to 580
1610 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1611  530 CALL xerrwd('SLSODE- At T(=R1) and step size H(=R2), the error',
1612  1 50, 204, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1613  CALL xerrwd(' test failed repeatedly or with ABS(H) = HMIN',
1614  1 50, 204, 0, 0, 0, 0, 2, tn, h)
1615  istate = -4
1616  go to 560
1617 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1618  540 CALL xerrwd('SLSODE- At T (=R1) and step size H (=R2), the ',
1619  1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1620  CALL xerrwd(' corrector convergence failed repeatedly ',
1621  1 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1622  CALL xerrwd(' or with ABS(H) = HMIN ',
1623  1 30, 205, 0, 0, 0, 0, 2, tn, h)
1624  istate = -5
1625 C Compute IMXER if relevant. -------------------------------------------
1626  560 big = 0.0e0
1627  imxer = 1
1628  DO 570 i = 1,n
1629  SIZE = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1630  IF (big .GE. size) go to 570
1631  big = SIZE
1632  imxer = i
1633  570 CONTINUE
1634  iwork(16) = imxer
1635 C Set Y vector, T, and optional outputs. -------------------------------
1636  580 DO 590 i = 1,n
1637  590 y(i) = rwork(i+lyh-1)
1638  t = tn
1639  rwork(11) = hu
1640  rwork(12) = h
1641  rwork(13) = tn
1642  iwork(11) = nst
1643  iwork(12) = nfe
1644  iwork(13) = nje
1645  iwork(14) = nqu
1646  iwork(15) = nq
1647  RETURN
1648 C-----------------------------------------------------------------------
1649 C Block I.
1650 C The following block handles all error returns due to illegal input
1651 C (ISTATE = -3), as detected before calling the core integrator.
1652 C First the error message routine is called. If the illegal input
1653 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1654 C-----------------------------------------------------------------------
1655  601 CALL xerrwd('SLSODE- ISTATE (=I1) illegal ',
1656  1 30, 1, 0, 1, istate, 0, 0, 0.0e0, 0.0e0)
1657  IF (istate .LT. 0) go to 800
1658  go to 700
1659  602 CALL xerrwd('SLSODE- ITASK (=I1) illegal ',
1660  1 30, 2, 0, 1, itask, 0, 0, 0.0e0, 0.0e0)
1661  go to 700
1662  603 CALL xerrwd(.GT.'SLSODE- ISTATE 1 but SLSODE not initialized ',
1663  1 50, 3, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1664  go to 700
1665  604 CALL xerrwd(.LT.'SLSODE- NEQ (=I1) 1 ',
1666  1 30, 4, 0, 1, neq(1), 0, 0, 0.0e0, 0.0e0)
1667  go to 700
1668  605 CALL xerrwd('SLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ',
1669  1 50, 5, 0, 2, n, neq(1), 0, 0.0e0, 0.0e0)
1670  go to 700
1671  606 CALL xerrwd('SLSODE- ITOL (=I1) illegal ',
1672  1 30, 6, 0, 1, itol, 0, 0, 0.0e0, 0.0e0)
1673  go to 700
1674  607 CALL xerrwd('SLSODE- IOPT (=I1) illegal ',
1675  1 30, 7, 0, 1, iopt, 0, 0, 0.0e0, 0.0e0)
1676  go to 700
1677  608 CALL xerrwd('SLSODE- MF (=I1) illegal ',
1678  1 30, 8, 0, 1, mf, 0, 0, 0.0e0, 0.0e0)
1679  go to 700
1680  609 CALL xerrwd(.LT..GE.'SLSODE- ML (=I1) illegal.. 0 or NEQ (=I2)',
1681  1 50, 9, 0, 2, ml, neq(1), 0, 0.0e0, 0.0e0)
1682  go to 700
1683  610 CALL xerrwd(.LT..GE.'SLSODE- MU (=I1) illegal.. 0 or NEQ (=I2)',
1684  1 50, 10, 0, 2, mu, neq(1), 0, 0.0e0, 0.0e0)
1685  go to 700
1686  611 CALL xerrwd(.LT.'SLSODE- MAXORD (=I1) 0 ',
1687  1 30, 11, 0, 1, maxord, 0, 0, 0.0e0, 0.0e0)
1688  go to 700
1689  612 CALL xerrwd(.LT.'SLSODE- MXSTEP (=I1) 0 ',
1690  1 30, 12, 0, 1, mxstep, 0, 0, 0.0e0, 0.0e0)
1691  go to 700
1692  613 CALL xerrwd(.LT.'SLSODE- MXHNIL (=I1) 0 ',
1693  1 30, 13, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1694  go to 700
1695  614 CALL xerrwd('SLSODE- TOUT (=R1) behind T (=R2) ',
1696  1 40, 14, 0, 0, 0, 0, 2, tout, t)
1697  CALL xerrwd(' Integration direction is given by H0 (=R1) ',
1698  1 50, 14, 0, 0, 0, 0, 1, h0, 0.0e0)
1699  go to 700
1700  615 CALL xerrwd(.LT.'SLSODE- HMAX (=R1) 0.0 ',
1701  1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0e0)
1702  go to 700
1703  616 CALL xerrwd(.LT.'SLSODE- HMIN (=R1) 0.0 ',
1704  1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0e0)
1705  go to 700
1706  617 CALL xerrwd(
1707  1 'SLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)',
1708  1 60, 17, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1709  go to 700
1710  618 CALL xerrwd(
1711  1 'SLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)',
1712  1 60, 18, 0, 2, leniw, liw, 0, 0.0e0, 0.0e0)
1713  go to 700
1714  619 CALL xerrwd(.LT.'SLSODE- RTOL(I1) is R1 0.0 ',
1715  1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0e0)
1716  go to 700
1717  620 CALL xerrwd(.LT.'SLSODE- ATOL(I1) is R1 0.0 ',
1718  1 40, 20, 0, 1, i, 0, 1, atoli, 0.0e0)
1719  go to 700
1720  621 ewti = rwork(lewt+i-1)
1721  CALL xerrwd(.LE.'SLSODE- EWT(I1) is R1 0.0 ',
1722  1 40, 21, 0, 1, i, 0, 1, ewti, 0.0e0)
1723  go to 700
1724  622 CALL xerrwd(
1725  1 'SLSODE- TOUT (=R1) too close to T(=R2) to start integration',
1726  1 60, 22, 0, 0, 0, 0, 2, tout, t)
1727  go to 700
1728  623 CALL xerrwd(
1729  1 'SLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ',
1730  1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1731  go to 700
1732  624 CALL xerrwd(
1733  1 'SLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ',
1734  1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1735  go to 700
1736  625 CALL xerrwd(
1737  1 'SLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ',
1738  1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1739  go to 700
1740  626 CALL xerrwd('SLSODE- At start of problem, too much accuracy ',
1741  1 50, 26, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1742  CALL xerrwd(
1743  1 ' requested for precision of machine.. See TOLSF (=R1) ',
1744  1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0e0)
1745  rwork(14) = tolsf
1746  go to 700
1747  627 CALL xerrwd('SLSODE- Trouble in SINTDY. ITASK = I1, TOUT = R1',
1748  1 50, 27, 0, 1, itask, 0, 1, tout, 0.0e0)
1749 C
1750  700 istate = -3
1751  RETURN
1752 C
1753  800 CALL xerrwd('SLSODE- Run aborted.. apparent infinite loop ',
1754  1 50, 303, 2, 0, 0, 0, 0, 0.0e0, 0.0e0)
1755  RETURN
1756 C----------------------- END OF SUBROUTINE SLSODE ----------------------
1757  END
ssolsy
subroutine ssolsy(WM, IWM, X, TEM)
Definition: ssolsy.f:1
sewset
subroutine sewset(N, ITOL, RTOL, ATOL, YCUR, EWT)
Definition: sewset.f:1
workspace_element::dimension
std::string dimension(void) const
Definition: workspace-element.h:74
xerrwd
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:3
slsode
subroutine slsode(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
Definition: slsode.f:2
sintdy
subroutine sintdy(T, K, YH, NYH, DKY, IFLAG)
Definition: sintdy.f:1
f
F77_RET_T const double const double * f
Definition: __pchip_deriv__.cc:38
size
size_t size(T const (&)[z])
Definition: help.cc:103
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
sprepj
subroutine sprepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
Definition: sprepj.f:1
sstode
subroutine sstode(NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
Definition: sstode.f:1
abs
T abs(T x)
Definition: pr-output.cc:3062
octave_value::sqrt
octave_value sqrt(void) const
Definition: ov.h:1200
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by   doxygen 1.8.8