| blob | 5a867955138d37bc611b492cc49222da6944b88a |
1 *DECK DLSODES
8 C-----------------------------------------------------------------------
9 C This is the 12 November 2003 version of
10 C DLSODES: Livermore Solver for Ordinary Differential Equations
11 C with general Sparse Jacobian matrix.
12 C
13 C This version is in double precision.
14 C
15 C DLSODES solves the initial value problem for stiff or nonstiff
16 C systems of first order ODEs,
17 C dy/dt = f(t,y) , or, in component form,
18 C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
19 C DLSODES is a variant of the DLSODE package, and is intended for
20 C problems in which the Jacobian matrix df/dy has an arbitrary
21 C sparse structure (when the problem is stiff).
22 C
23 C Authors: Alan C. Hindmarsh
24 C Center for Applied Scientific Computing, L-561
25 C Lawrence Livermore National Laboratory
26 C Livermore, CA 94551
27 C and
28 C Andrew H. Sherman
29 C J. S. Nolen and Associates
30 C Houston, TX 77084
31 C-----------------------------------------------------------------------
32 C References:
33 C 1. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE
34 C Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.),
35 C North-Holland, Amsterdam, 1983, pp. 55-64.
36 C
37 C 2. S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman,
38 C Yale Sparse Matrix Package: I. The Symmetric Codes,
39 C Int. J. Num. Meth. Eng., 18 (1982), pp. 1145-1151.
40 C
41 C 3. S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman,
42 C Yale Sparse Matrix Package: II. The Nonsymmetric Codes,
43 C Research Report No. 114, Dept. of Computer Sciences, Yale
44 C University, 1977.
45 C-----------------------------------------------------------------------
46 C Summary of Usage.
47 C
48 C Communication between the user and the DLSODES package, for normal
49 C situations, is summarized here. This summary describes only a subset
50 C of the full set of options available. See the full description for
51 C details, including optional communication, nonstandard options,
52 C and instructions for special situations. See also the example
53 C problem (with program and output) following this summary.
54 C
55 C A. First provide a subroutine of the form:
56 C SUBROUTINE F (NEQ, T, Y, YDOT)
57 C DOUBLE PRECISION T, Y(*), YDOT(*)
58 C which supplies the vector function f by loading YDOT(i) with f(i).
59 C
60 C B. Next determine (or guess) whether or not the problem is stiff.
61 C Stiffness occurs when the Jacobian matrix df/dy has an eigenvalue
62 C whose real part is negative and large in magnitude, compared to the
63 C reciprocal of the t span of interest. If the problem is nonstiff,
64 C use a method flag MF = 10. If it is stiff, there are two standard
65 C choices for the method flag, MF = 121 and MF = 222. In both cases,
66 C DLSODES requires the Jacobian matrix in some form, and it treats this
67 C matrix in general sparse form, with sparsity structure determined
68 C internally. (For options where the user supplies the sparsity
69 C structure, see the full description of MF below.)
70 C
71 C C. If the problem is stiff, you are encouraged to supply the Jacobian
72 C directly (MF = 121), but if this is not feasible, DLSODES will
73 C compute it internally by difference quotients (MF = 222).
74 C If you are supplying the Jacobian, provide a subroutine of the form:
75 C SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ)
76 C DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*)
77 C Here NEQ, T, Y, and J are input arguments, and the JAC routine is to
78 C load the array PDJ (of length NEQ) with the J-th column of df/dy.
79 C I.e., load PDJ(i) with df(i)/dy(J) for all relevant values of i.
80 C The arguments IAN and JAN should be ignored for normal situations.
81 C DLSODES will call the JAC routine with J = 1,2,...,NEQ.
82 C Only nonzero elements need be loaded. Usually, a crude approximation
83 C to df/dy, possibly with fewer nonzero elements, will suffice.
84 C
85 C D. Write a main program which calls Subroutine DLSODES once for
86 C each point at which answers are desired. This should also provide
87 C for possible use of logical unit 6 for output of error messages by
88 C DLSODES. On the first call to DLSODES, supply arguments as follows:
89 C F = name of subroutine for right-hand side vector f.
90 C This name must be declared External in calling program.
91 C NEQ = number of first order ODEs.
92 C Y = array of initial values, of length NEQ.
93 C T = the initial value of the independent variable t.
94 C TOUT = first point where output is desired (.ne. T).
95 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
96 C RTOL = relative tolerance parameter (scalar).
97 C ATOL = absolute tolerance parameter (scalar or array).
98 C The estimated local error in Y(i) will be controlled so as
99 C to be roughly less (in magnitude) than
100 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
101 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
102 C Thus the local error test passes if, in each component,
103 C either the absolute error is less than ATOL (or ATOL(i)),
104 C or the relative error is less than RTOL.
105 C Use RTOL = 0.0 for pure absolute error control, and
106 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
107 C control. Caution: actual (global) errors may exceed these
108 C local tolerances, so choose them conservatively.
109 C ITASK = 1 for normal computation of output values of Y at t = TOUT.
110 C ISTATE = integer flag (input and output). Set ISTATE = 1.
111 C IOPT = 0 to indicate no optional inputs used.
112 C RWORK = real work array of length at least:
113 C 20 + 16*NEQ for MF = 10,
114 C 20 + (2 + 1./LENRAT)*NNZ + (11 + 9./LENRAT)*NEQ
115 C for MF = 121 or 222,
116 C where:
117 C NNZ = the number of nonzero elements in the sparse
118 C Jacobian (if this is unknown, use an estimate), and
119 C LENRAT = the real to integer wordlength ratio (usually 1 in
120 C single precision and 2 in double precision).
121 C In any case, the required size of RWORK cannot generally
122 C be predicted in advance if MF = 121 or 222, and the value
123 C above is a rough estimate of a crude lower bound. Some
124 C experimentation with this size may be necessary.
125 C (When known, the correct required length is an optional
126 C output, available in IWORK(17).)
127 C LRW = declared length of RWORK (in user dimension).
128 C IWORK = integer work array of length at least 30.
129 C LIW = declared length of IWORK (in user dimension).
130 C JAC = name of subroutine for Jacobian matrix (MF = 121).
131 C If used, this name must be declared External in calling
132 C program. If not used, pass a dummy name.
133 C MF = method flag. Standard values are:
134 C 10 for nonstiff (Adams) method, no Jacobian used
135 C 121 for stiff (BDF) method, user-supplied sparse Jacobian
136 C 222 for stiff method, internally generated sparse Jacobian
137 C Note that the main program must declare arrays Y, RWORK, IWORK,
138 C and possibly ATOL.
139 C
140 C E. The output from the first call (or any call) is:
141 C Y = array of computed values of y(t) vector.
142 C T = corresponding value of independent variable (normally TOUT).
143 C ISTATE = 2 if DLSODES was successful, negative otherwise.
144 C -1 means excess work done on this call (perhaps wrong MF).
145 C -2 means excess accuracy requested (tolerances too small).
146 C -3 means illegal input detected (see printed message).
147 C -4 means repeated error test failures (check all inputs).
148 C -5 means repeated convergence failures (perhaps bad Jacobian
149 C supplied or wrong choice of MF or tolerances).
150 C -6 means error weight became zero during problem. (Solution
151 C component i vanished, and ATOL or ATOL(i) = 0.)
152 C -7 means a fatal error return flag came from sparse solver
153 C CDRV by way of DPRJS or DSOLSS. Should never happen.
154 C A return with ISTATE = -1, -4, or -5 may result from using
155 C an inappropriate sparsity structure, one that is quite
156 C different from the initial structure. Consider calling
157 C DLSODES again with ISTATE = 3 to force the structure to be
158 C reevaluated. See the full description of ISTATE below.
159 C
160 C F. To continue the integration after a successful return, simply
161 C reset TOUT and call DLSODES again. No other parameters need be reset.
162 C
163 C-----------------------------------------------------------------------
164 C Example Problem.
165 C
166 C The following is a simple example problem, with the coding
167 C needed for its solution by DLSODES. The problem is from chemical
168 C kinetics, and consists of the following 12 rate equations:
169 C dy1/dt = -rk1*y1
170 C dy2/dt = rk1*y1 + rk11*rk14*y4 + rk19*rk14*y5
171 C - rk3*y2*y3 - rk15*y2*y12 - rk2*y2
172 C dy3/dt = rk2*y2 - rk5*y3 - rk3*y2*y3 - rk7*y10*y3
173 C + rk11*rk14*y4 + rk12*rk14*y6
174 C dy4/dt = rk3*y2*y3 - rk11*rk14*y4 - rk4*y4
175 C dy5/dt = rk15*y2*y12 - rk19*rk14*y5 - rk16*y5
176 C dy6/dt = rk7*y10*y3 - rk12*rk14*y6 - rk8*y6
177 C dy7/dt = rk17*y10*y12 - rk20*rk14*y7 - rk18*y7
178 C dy8/dt = rk9*y10 - rk13*rk14*y8 - rk10*y8
179 C dy9/dt = rk4*y4 + rk16*y5 + rk8*y6 + rk18*y7
180 C dy10/dt = rk5*y3 + rk12*rk14*y6 + rk20*rk14*y7
181 C + rk13*rk14*y8 - rk7*y10*y3 - rk17*y10*y12
182 C - rk6*y10 - rk9*y10
183 C dy11/dt = rk10*y8
184 C dy12/dt = rk6*y10 + rk19*rk14*y5 + rk20*rk14*y7
185 C - rk15*y2*y12 - rk17*y10*y12
186 C
187 C with rk1 = rk5 = 0.1, rk4 = rk8 = rk16 = rk18 = 2.5,
188 C rk10 = 5.0, rk2 = rk6 = 10.0, rk14 = 30.0,
189 C rk3 = rk7 = rk9 = rk11 = rk12 = rk13 = rk19 = rk20 = 50.0,
190 C rk15 = rk17 = 100.0.
191 C
192 C The t interval is from 0 to 1000, and the initial conditions
193 C are y1 = 1, y2 = y3 = ... = y12 = 0. The problem is stiff.
194 C
195 C The following coding solves this problem with DLSODES, using MF = 121
196 C and printing results at t = .1, 1., 10., 100., 1000. It uses
197 C ITOL = 1 and mixed relative/absolute tolerance controls.
198 C During the run and at the end, statistical quantities of interest
199 C are printed (see optional outputs in the full description below).
200 C
201 C EXTERNAL FEX, JEX
202 C DOUBLE PRECISION ATOL, RTOL, RWORK, T, TOUT, Y
203 C DIMENSION Y(12), RWORK(500), IWORK(30)
204 C DATA LRW/500/, LIW/30/
205 C NEQ = 12
206 C DO 10 I = 1,NEQ
207 C 10 Y(I) = 0.0D0
208 C Y(1) = 1.0D0
209 C T = 0.0D0
210 C TOUT = 0.1D0
211 C ITOL = 1
212 C RTOL = 1.0D-4
213 C ATOL = 1.0D-6
214 C ITASK = 1
215 C ISTATE = 1
216 C IOPT = 0
217 C MF = 121
218 C DO 40 IOUT = 1,5
219 C CALL DLSODES (FEX, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL,
220 C 1 ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
221 C WRITE(6,30)T,IWORK(11),RWORK(11),(Y(I),I=1,NEQ)
222 C 30 FORMAT(//' At t =',D11.3,4X,
223 C 1 ' No. steps =',I5,4X,' Last step =',D11.3/
224 C 2 ' Y array = ',4D14.5/13X,4D14.5/13X,4D14.5)
225 C IF (ISTATE .LT. 0) GO TO 80
226 C TOUT = TOUT*10.0D0
227 C 40 CONTINUE
228 C LENRW = IWORK(17)
229 C LENIW = IWORK(18)
230 C NST = IWORK(11)
231 C NFE = IWORK(12)
232 C NJE = IWORK(13)
233 C NLU = IWORK(21)
234 C NNZ = IWORK(19)
235 C NNZLU = IWORK(25) + IWORK(26) + NEQ
236 C WRITE (6,70) LENRW,LENIW,NST,NFE,NJE,NLU,NNZ,NNZLU
237 C 70 FORMAT(//' Required RWORK size =',I4,' IWORK size =',I4/
238 C 1 ' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4,
239 C 2 ' No. LU-s =',I4/' No. of nonzeros in J =',I5,
240 C 3 ' No. of nonzeros in LU =',I5)
241 C STOP
242 C 80 WRITE(6,90)ISTATE
243 C 90 FORMAT(///' Error halt.. ISTATE =',I3)
244 C STOP
245 C END
246 C
247 C SUBROUTINE FEX (NEQ, T, Y, YDOT)
248 C DOUBLE PRECISION T, Y, YDOT
249 C DOUBLE PRECISION RK1, RK2, RK3, RK4, RK5, RK6, RK7, RK8, RK9,
250 C 1 RK10, RK11, RK12, RK13, RK14, RK15, RK16, RK17
251 C DIMENSION Y(12), YDOT(12)
252 C DATA RK1/0.1D0/, RK2/10.0D0/, RK3/50.0D0/, RK4/2.5D0/, RK5/0.1D0/,
253 C 1 RK6/10.0D0/, RK7/50.0D0/, RK8/2.5D0/, RK9/50.0D0/, RK10/5.0D0/,
254 C 2 RK11/50.0D0/, RK12/50.0D0/, RK13/50.0D0/, RK14/30.0D0/,
255 C 3 RK15/100.0D0/, RK16/2.5D0/, RK17/100.0D0/, RK18/2.5D0/,
256 C 4 RK19/50.0D0/, RK20/50.0D0/
257 C YDOT(1) = -RK1*Y(1)
258 C YDOT(2) = RK1*Y(1) + RK11*RK14*Y(4) + RK19*RK14*Y(5)
259 C 1 - RK3*Y(2)*Y(3) - RK15*Y(2)*Y(12) - RK2*Y(2)
260 C YDOT(3) = RK2*Y(2) - RK5*Y(3) - RK3*Y(2)*Y(3) - RK7*Y(10)*Y(3)
261 C 1 + RK11*RK14*Y(4) + RK12*RK14*Y(6)
262 C YDOT(4) = RK3*Y(2)*Y(3) - RK11*RK14*Y(4) - RK4*Y(4)
263 C YDOT(5) = RK15*Y(2)*Y(12) - RK19*RK14*Y(5) - RK16*Y(5)
264 C YDOT(6) = RK7*Y(10)*Y(3) - RK12*RK14*Y(6) - RK8*Y(6)
265 C YDOT(7) = RK17*Y(10)*Y(12) - RK20*RK14*Y(7) - RK18*Y(7)
266 C YDOT(8) = RK9*Y(10) - RK13*RK14*Y(8) - RK10*Y(8)
267 C YDOT(9) = RK4*Y(4) + RK16*Y(5) + RK8*Y(6) + RK18*Y(7)
268 C YDOT(10) = RK5*Y(3) + RK12*RK14*Y(6) + RK20*RK14*Y(7)
269 C 1 + RK13*RK14*Y(8) - RK7*Y(10)*Y(3) - RK17*Y(10)*Y(12)
270 C 2 - RK6*Y(10) - RK9*Y(10)
271 C YDOT(11) = RK10*Y(8)
272 C YDOT(12) = RK6*Y(10) + RK19*RK14*Y(5) + RK20*RK14*Y(7)
273 C 1 - RK15*Y(2)*Y(12) - RK17*Y(10)*Y(12)
274 C RETURN
275 C END
276 C
277 C SUBROUTINE JEX (NEQ, T, Y, J, IA, JA, PDJ)
278 C DOUBLE PRECISION T, Y, PDJ
279 C DOUBLE PRECISION RK1, RK2, RK3, RK4, RK5, RK6, RK7, RK8, RK9,
280 C 1 RK10, RK11, RK12, RK13, RK14, RK15, RK16, RK17
281 C DIMENSION Y(12), IA(*), JA(*), PDJ(12)
282 C DATA RK1/0.1D0/, RK2/10.0D0/, RK3/50.0D0/, RK4/2.5D0/, RK5/0.1D0/,
283 C 1 RK6/10.0D0/, RK7/50.0D0/, RK8/2.5D0/, RK9/50.0D0/, RK10/5.0D0/,
284 C 2 RK11/50.0D0/, RK12/50.0D0/, RK13/50.0D0/, RK14/30.0D0/,
285 C 3 RK15/100.0D0/, RK16/2.5D0/, RK17/100.0D0/, RK18/2.5D0/,
286 C 4 RK19/50.0D0/, RK20/50.0D0/
287 C GO TO (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), J
288 C 1 PDJ(1) = -RK1
289 C PDJ(2) = RK1
290 C RETURN
291 C 2 PDJ(2) = -RK3*Y(3) - RK15*Y(12) - RK2
292 C PDJ(3) = RK2 - RK3*Y(3)
293 C PDJ(4) = RK3*Y(3)
294 C PDJ(5) = RK15*Y(12)
295 C PDJ(12) = -RK15*Y(12)
296 C RETURN
297 C 3 PDJ(2) = -RK3*Y(2)
298 C PDJ(3) = -RK5 - RK3*Y(2) - RK7*Y(10)
299 C PDJ(4) = RK3*Y(2)
300 C PDJ(6) = RK7*Y(10)
301 C PDJ(10) = RK5 - RK7*Y(10)
302 C RETURN
303 C 4 PDJ(2) = RK11*RK14
304 C PDJ(3) = RK11*RK14
305 C PDJ(4) = -RK11*RK14 - RK4
306 C PDJ(9) = RK4
307 C RETURN
308 C 5 PDJ(2) = RK19*RK14
309 C PDJ(5) = -RK19*RK14 - RK16
310 C PDJ(9) = RK16
311 C PDJ(12) = RK19*RK14
312 C RETURN
313 C 6 PDJ(3) = RK12*RK14
314 C PDJ(6) = -RK12*RK14 - RK8
315 C PDJ(9) = RK8
316 C PDJ(10) = RK12*RK14
317 C RETURN
318 C 7 PDJ(7) = -RK20*RK14 - RK18
319 C PDJ(9) = RK18
320 C PDJ(10) = RK20*RK14
321 C PDJ(12) = RK20*RK14
322 C RETURN
323 C 8 PDJ(8) = -RK13*RK14 - RK10
324 C PDJ(10) = RK13*RK14
325 C PDJ(11) = RK10
326 C 9 RETURN
327 C 10 PDJ(3) = -RK7*Y(3)
328 C PDJ(6) = RK7*Y(3)
329 C PDJ(7) = RK17*Y(12)
330 C PDJ(8) = RK9
331 C PDJ(10) = -RK7*Y(3) - RK17*Y(12) - RK6 - RK9
332 C PDJ(12) = RK6 - RK17*Y(12)
333 C 11 RETURN
334 C 12 PDJ(2) = -RK15*Y(2)
335 C PDJ(5) = RK15*Y(2)
336 C PDJ(7) = RK17*Y(10)
337 C PDJ(10) = -RK17*Y(10)
338 C PDJ(12) = -RK15*Y(2) - RK17*Y(10)
339 C RETURN
340 C END
341 C
342 C The output of this program (on a Cray-1 in single precision)
343 C is as follows:
344 C
345 C
346 C At t = 1.000e-01 No. steps = 12 Last step = 1.515e-02
347 C Y array = 9.90050e-01 6.28228e-03 3.65313e-03 7.51934e-07
348 C 1.12167e-09 1.18458e-09 1.77291e-12 3.26476e-07
349 C 5.46720e-08 9.99500e-06 4.48483e-08 2.76398e-06
350 C
351 C
352 C At t = 1.000e+00 No. steps = 33 Last step = 7.880e-02
353 C Y array = 9.04837e-01 9.13105e-03 8.20622e-02 2.49177e-05
354 C 1.85055e-06 1.96797e-06 1.46157e-07 2.39557e-05
355 C 3.26306e-05 7.21621e-04 5.06433e-05 3.05010e-03
356 C
357 C
358 C At t = 1.000e+01 No. steps = 48 Last step = 1.239e+00
359 C Y array = 3.67876e-01 3.68958e-03 3.65133e-01 4.48325e-05
360 C 6.10798e-05 4.33148e-05 5.90211e-05 1.18449e-04
361 C 3.15235e-03 3.56531e-03 4.15520e-03 2.48741e-01
362 C
363 C
364 C At t = 1.000e+02 No. steps = 91 Last step = 3.764e+00
365 C Y array = 4.44981e-05 4.42666e-07 4.47273e-04 -3.53257e-11
366 C 2.81577e-08 -9.67741e-11 2.77615e-07 1.45322e-07
367 C 1.56230e-02 4.37394e-06 1.60104e-02 9.52246e-01
368 C
369 C
370 C At t = 1.000e+03 No. steps = 111 Last step = 4.156e+02
371 C Y array = -2.65492e-13 2.60539e-14 -8.59563e-12 6.29355e-14
372 C -1.78066e-13 5.71471e-13 -1.47561e-12 4.58078e-15
373 C 1.56314e-02 1.37878e-13 1.60184e-02 9.52719e-01
374 C
375 C
376 C Required RWORK size = 442 IWORK size = 30
377 C No. steps = 111 No. f-s = 142 No. J-s = 2 No. LU-s = 20
378 C No. of nonzeros in J = 44 No. of nonzeros in LU = 50
379 C
380 C-----------------------------------------------------------------------
381 C Full Description of User Interface to DLSODES.
382 C
383 C The user interface to DLSODES consists of the following parts.
384 C
385 C 1. The call sequence to Subroutine DLSODES, which is a driver
386 C routine for the solver. This includes descriptions of both
387 C the call sequence arguments and of user-supplied routines.
388 C Following these descriptions is a description of
389 C optional inputs available through the call sequence, and then
390 C a description of optional outputs (in the work arrays).
391 C
392 C 2. Descriptions of other routines in the DLSODES package that may be
393 C (optionally) called by the user. These provide the ability to
394 C alter error message handling, save and restore the internal
395 C Common, and obtain specified derivatives of the solution y(t).
396 C
397 C 3. Descriptions of Common blocks to be declared in overlay
398 C or similar environments, or to be saved when doing an interrupt
399 C of the problem and continued solution later.
400 C
401 C 4. Description of two routines in the DLSODES package, either of
402 C which the user may replace with his/her own version, if desired.
403 C These relate to the measurement of errors.
404 C
405 C-----------------------------------------------------------------------
406 C Part 1. Call Sequence.
407 C
408 C The call sequence parameters used for input only are
409 C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
410 C and those used for both input and output are
411 C Y, T, ISTATE.
412 C The work arrays RWORK and IWORK are also used for conditional and
413 C optional inputs and optional outputs. (The term output here refers
414 C to the return from Subroutine DLSODES to the user's calling program.)
415 C
416 C The legality of input parameters will be thoroughly checked on the
417 C initial call for the problem, but not checked thereafter unless a
418 C change in input parameters is flagged by ISTATE = 3 on input.
419 C
420 C The descriptions of the call arguments are as follows.
421 C
422 C F = the name of the user-supplied subroutine defining the
423 C ODE system. The system must be put in the first-order
424 C form dy/dt = f(t,y), where f is a vector-valued function
425 C of the scalar t and the vector y. Subroutine F is to
426 C compute the function f. It is to have the form
427 C SUBROUTINE F (NEQ, T, Y, YDOT)
428 C DOUBLE PRECISION T, Y(*), YDOT(*)
429 C where NEQ, T, and Y are input, and the array YDOT = f(t,y)
430 C is output. Y and YDOT are arrays of length NEQ.
431 C Subroutine F should not alter y(1),...,y(NEQ).
432 C F must be declared External in the calling program.
433 C
434 C Subroutine F may access user-defined quantities in
435 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
436 C (dimensioned in F) and/or Y has length exceeding NEQ(1).
437 C See the descriptions of NEQ and Y below.
438 C
439 C If quantities computed in the F routine are needed
440 C externally to DLSODES, an extra call to F should be made
441 C for this purpose, for consistent and accurate results.
442 C If only the derivative dy/dt is needed, use DINTDY instead.
443 C
444 C NEQ = the size of the ODE system (number of first order
445 C ordinary differential equations). Used only for input.
446 C NEQ may be decreased, but not increased, during the problem.
447 C If NEQ is decreased (with ISTATE = 3 on input), the
448 C remaining components of Y should be left undisturbed, if
449 C these are to be accessed in F and/or JAC.
450 C
451 C Normally, NEQ is a scalar, and it is generally referred to
452 C as a scalar in this user interface description. However,
453 C NEQ may be an array, with NEQ(1) set to the system size.
454 C (The DLSODES package accesses only NEQ(1).) In either case,
455 C this parameter is passed as the NEQ argument in all calls
456 C to F and JAC. Hence, if it is an array, locations
457 C NEQ(2),... may be used to store other integer data and pass
458 C it to F and/or JAC. Subroutines F and/or JAC must include
459 C NEQ in a Dimension statement in that case.
460 C
461 C Y = a real array for the vector of dependent variables, of
462 C length NEQ or more. Used for both input and output on the
463 C first call (ISTATE = 1), and only for output on other calls.
464 C on the first call, Y must contain the vector of initial
465 C values. On output, Y contains the computed solution vector,
466 C evaluated at T. If desired, the Y array may be used
467 C for other purposes between calls to the solver.
468 C
469 C This array is passed as the Y argument in all calls to
470 C F and JAC. Hence its length may exceed NEQ, and locations
471 C Y(NEQ+1),... may be used to store other real data and
472 C pass it to F and/or JAC. (The DLSODES package accesses only
473 C Y(1),...,Y(NEQ).)
474 C
475 C T = the independent variable. On input, T is used only on the
476 C first call, as the initial point of the integration.
477 C on output, after each call, T is the value at which a
478 C computed solution Y is evaluated (usually the same as TOUT).
479 C On an error return, T is the farthest point reached.
480 C
481 C TOUT = the next value of t at which a computed solution is desired.
482 C Used only for input.
483 C
484 C When starting the problem (ISTATE = 1), TOUT may be equal
485 C to T for one call, then should .ne. T for the next call.
486 C For the initial T, an input value of TOUT .ne. T is used
487 C in order to determine the direction of the integration
488 C (i.e. the algebraic sign of the step sizes) and the rough
489 C scale of the problem. Integration in either direction
490 C (forward or backward in t) is permitted.
491 C
492 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
493 C the first call (i.e. the first call with TOUT .ne. T).
494 C Otherwise, TOUT is required on every call.
495 C
496 C If ITASK = 1, 3, or 4, the values of TOUT need not be
497 C monotone, but a value of TOUT which backs up is limited
498 C to the current internal T interval, whose endpoints are
499 C TCUR - HU and TCUR (see optional outputs, below, for
500 C TCUR and HU).
501 C
502 C ITOL = an indicator for the type of error control. See
503 C description below under ATOL. Used only for input.
504 C
505 C RTOL = a relative error tolerance parameter, either a scalar or
506 C an array of length NEQ. See description below under ATOL.
507 C Input only.
508 C
509 C ATOL = an absolute error tolerance parameter, either a scalar or
510 C an array of length NEQ. Input only.
511 C
512 C The input parameters ITOL, RTOL, and ATOL determine
513 C the error control performed by the solver. The solver will
514 C control the vector E = (E(i)) of estimated local errors
515 C in y, according to an inequality of the form
516 C RMS-norm of ( E(i)/EWT(i) ) .le. 1,
517 C where EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
518 C and the RMS-norm (root-mean-square norm) here is
519 C RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
520 C is a vector of weights which must always be positive, and
521 C the values of RTOL and ATOL should all be non-negative.
522 C The following table gives the types (scalar/array) of
523 C RTOL and ATOL, and the corresponding form of EWT(i).
524 C
525 C ITOL RTOL ATOL EWT(i)
526 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
527 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
528 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
529 C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i)
530 C
531 C When either of these parameters is a scalar, it need not
532 C be dimensioned in the user's calling program.
533 C
534 C If none of the above choices (with ITOL, RTOL, and ATOL
535 C fixed throughout the problem) is suitable, more general
536 C error controls can be obtained by substituting
537 C user-supplied routines for the setting of EWT and/or for
538 C the norm calculation. See Part 4 below.
539 C
540 C If global errors are to be estimated by making a repeated
541 C run on the same problem with smaller tolerances, then all
542 C components of RTOL and ATOL (i.e. of EWT) should be scaled
543 C down uniformly.
544 C
545 C ITASK = an index specifying the task to be performed.
546 C Input only. ITASK has the following values and meanings.
547 C 1 means normal computation of output values of y(t) at
548 C t = TOUT (by overshooting and interpolating).
549 C 2 means take one step only and return.
550 C 3 means stop at the first internal mesh point at or
551 C beyond t = TOUT and return.
552 C 4 means normal computation of output values of y(t) at
553 C t = TOUT but without overshooting t = TCRIT.
554 C TCRIT must be input as RWORK(1). TCRIT may be equal to
555 C or beyond TOUT, but not behind it in the direction of
556 C integration. This option is useful if the problem
557 C has a singularity at or beyond t = TCRIT.
558 C 5 means take one step, without passing TCRIT, and return.
559 C TCRIT must be input as RWORK(1).
560 C
561 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
562 C (within roundoff), it will return T = TCRIT (exactly) to
563 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
564 C in which case answers at t = TOUT are returned first).
565 C
566 C ISTATE = an index used for input and output to specify the
567 C the state of the calculation.
568 C
569 C On input, the values of ISTATE are as follows.
570 C 1 means this is the first call for the problem
571 C (initializations will be done). See note below.
572 C 2 means this is not the first call, and the calculation
573 C is to continue normally, with no change in any input
574 C parameters except possibly TOUT and ITASK.
575 C (If ITOL, RTOL, and/or ATOL are changed between calls
576 C with ISTATE = 2, the new values will be used but not
577 C tested for legality.)
578 C 3 means this is not the first call, and the
579 C calculation is to continue normally, but with
580 C a change in input parameters other than
581 C TOUT and ITASK. Changes are allowed in
582 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF,
583 C the conditional inputs IA and JA,
584 C and any of the optional inputs except H0.
585 C In particular, if MITER = 1 or 2, a call with ISTATE = 3
586 C will cause the sparsity structure of the problem to be
587 C recomputed (or reread from IA and JA if MOSS = 0).
588 C Note: a preliminary call with TOUT = T is not counted
589 C as a first call here, as no initialization or checking of
590 C input is done. (Such a call is sometimes useful for the
591 C purpose of outputting the initial conditions.)
592 C Thus the first call for which TOUT .ne. T requires
593 C ISTATE = 1 on input.
594 C
595 C On output, ISTATE has the following values and meanings.
596 C 1 means nothing was done; TOUT = T and ISTATE = 1 on input.
597 C 2 means the integration was performed successfully.
598 C -1 means an excessive amount of work (more than MXSTEP
599 C steps) was done on this call, before completing the
600 C requested task, but the integration was otherwise
601 C successful as far as T. (MXSTEP is an optional input
602 C and is normally 500.) To continue, the user may
603 C simply reset ISTATE to a value .gt. 1 and call again
604 C (the excess work step counter will be reset to 0).
605 C In addition, the user may increase MXSTEP to avoid
606 C this error return (see below on optional inputs).
607 C -2 means too much accuracy was requested for the precision
608 C of the machine being used. This was detected before
609 C completing the requested task, but the integration
610 C was successful as far as T. To continue, the tolerance
611 C parameters must be reset, and ISTATE must be set
612 C to 3. The optional output TOLSF may be used for this
613 C purpose. (Note: If this condition is detected before
614 C taking any steps, then an illegal input return
615 C (ISTATE = -3) occurs instead.)
616 C -3 means illegal input was detected, before taking any
617 C integration steps. See written message for details.
618 C Note: If the solver detects an infinite loop of calls
619 C to the solver with illegal input, it will cause
620 C the run to stop.
621 C -4 means there were repeated error test failures on
622 C one attempted step, before completing the requested
623 C task, but the integration was successful as far as T.
624 C The problem may have a singularity, or the input
625 C may be inappropriate.
626 C -5 means there were repeated convergence test failures on
627 C one attempted step, before completing the requested
628 C task, but the integration was successful as far as T.
629 C This may be caused by an inaccurate Jacobian matrix,
630 C if one is being used.
631 C -6 means EWT(i) became zero for some i during the
632 C integration. Pure relative error control (ATOL(i)=0.0)
633 C was requested on a variable which has now vanished.
634 C The integration was successful as far as T.
635 C -7 means a fatal error return flag came from the sparse
636 C solver CDRV by way of DPRJS or DSOLSS (numerical
637 C factorization or backsolve). This should never happen.
638 C The integration was successful as far as T.
639 C
640 C Note: an error return with ISTATE = -1, -4, or -5 and with
641 C MITER = 1 or 2 may mean that the sparsity structure of the
642 C problem has changed significantly since it was last
643 C determined (or input). In that case, one can attempt to
644 C complete the integration by setting ISTATE = 3 on the next
645 C call, so that a new structure determination is done.
646 C
647 C Note: since the normal output value of ISTATE is 2,
648 C it does not need to be reset for normal continuation.
649 C Also, since a negative input value of ISTATE will be
650 C regarded as illegal, a negative output value requires the
651 C user to change it, and possibly other inputs, before
652 C calling the solver again.
653 C
654 C IOPT = an integer flag to specify whether or not any optional
655 C inputs are being used on this call. Input only.
656 C The optional inputs are listed separately below.
657 C IOPT = 0 means no optional inputs are being used.
658 C Default values will be used in all cases.
659 C IOPT = 1 means one or more optional inputs are being used.
660 C
661 C RWORK = a work array used for a mixture of real (double precision)
662 C and integer work space.
663 C The length of RWORK (in real words) must be at least
664 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM where
665 C NYH = the initial value of NEQ,
666 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
667 C smaller value is given as an optional input),
668 C LWM = 0 if MITER = 0,
669 C LWM = 2*NNZ + 2*NEQ + (NNZ+9*NEQ)/LENRAT if MITER = 1,
670 C LWM = 2*NNZ + 2*NEQ + (NNZ+10*NEQ)/LENRAT if MITER = 2,
671 C LWM = NEQ + 2 if MITER = 3.
672 C In the above formulas,
673 C NNZ = number of nonzero elements in the Jacobian matrix.
674 C LENRAT = the real to integer wordlength ratio (usually 1 in
675 C single precision and 2 in double precision).
676 C (See the MF description for METH and MITER.)
677 C Thus if MAXORD has its default value and NEQ is constant,
678 C the minimum length of RWORK is:
679 C 20 + 16*NEQ for MF = 10,
680 C 20 + 16*NEQ + LWM for MF = 11, 111, 211, 12, 112, 212,
681 C 22 + 17*NEQ for MF = 13,
682 C 20 + 9*NEQ for MF = 20,
683 C 20 + 9*NEQ + LWM for MF = 21, 121, 221, 22, 122, 222,
684 C 22 + 10*NEQ for MF = 23.
685 C If MITER = 1 or 2, the above formula for LWM is only a
686 C crude lower bound. The required length of RWORK cannot
687 C be readily predicted in general, as it depends on the
688 C sparsity structure of the problem. Some experimentation
689 C may be necessary.
690 C
691 C The first 20 words of RWORK are reserved for conditional
692 C and optional inputs and optional outputs.
693 C
694 C The following word in RWORK is a conditional input:
695 C RWORK(1) = TCRIT = critical value of t which the solver
696 C is not to overshoot. Required if ITASK is
697 C 4 or 5, and ignored otherwise. (See ITASK.)
698 C
699 C LRW = the length of the array RWORK, as declared by the user.
700 C (This will be checked by the solver.)
701 C
702 C IWORK = an integer work array. The length of IWORK must be at least
703 C 31 + NEQ + NNZ if MOSS = 0 and MITER = 1 or 2, or
704 C 30 otherwise.
705 C (NNZ is the number of nonzero elements in df/dy.)
706 C
707 C In DLSODES, IWORK is used only for conditional and
708 C optional inputs and optional outputs.
709 C
710 C The following two blocks of words in IWORK are conditional
711 C inputs, required if MOSS = 0 and MITER = 1 or 2, but not
712 C otherwise (see the description of MF for MOSS).
713 C IWORK(30+j) = IA(j) (j=1,...,NEQ+1)
714 C IWORK(31+NEQ+k) = JA(k) (k=1,...,NNZ)
715 C The two arrays IA and JA describe the sparsity structure
716 C to be assumed for the Jacobian matrix. JA contains the row
717 C indices where nonzero elements occur, reading in columnwise
718 C order, and IA contains the starting locations in JA of the
719 C descriptions of columns 1,...,NEQ, in that order, with
720 C IA(1) = 1. Thus, for each column index j = 1,...,NEQ, the
721 C values of the row index i in column j where a nonzero
722 C element may occur are given by
723 C i = JA(k), where IA(j) .le. k .lt. IA(j+1).
724 C If NNZ is the total number of nonzero locations assumed,
725 C then the length of the JA array is NNZ, and IA(NEQ+1) must
726 C be NNZ + 1. Duplicate entries are not allowed.
727 C
728 C LIW = the length of the array IWORK, as declared by the user.
729 C (This will be checked by the solver.)
730 C
731 C Note: The work arrays must not be altered between calls to DLSODES
732 C for the same problem, except possibly for the conditional and
733 C optional inputs, and except for the last 3*NEQ words of RWORK.
734 C The latter space is used for internal scratch space, and so is
735 C available for use by the user outside DLSODES between calls, if
736 C desired (but not for use by F or JAC).
737 C
738 C JAC = name of user-supplied routine (MITER = 1 or MOSS = 1) to
739 C compute the Jacobian matrix, df/dy, as a function of
740 C the scalar t and the vector y. It is to have the form
741 C SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ)
742 C DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*)
743 C where NEQ, T, Y, J, IAN, and JAN are input, and the array
744 C PDJ, of length NEQ, is to be loaded with column J
745 C of the Jacobian on output. Thus df(i)/dy(J) is to be
746 C loaded into PDJ(i) for all relevant values of i.
747 C Here T and Y have the same meaning as in Subroutine F,
748 C and J is a column index (1 to NEQ). IAN and JAN are
749 C undefined in calls to JAC for structure determination
750 C (MOSS = 1). otherwise, IAN and JAN are structure
751 C descriptors, as defined under optional outputs below, and
752 C so can be used to determine the relevant row indices i, if
753 C desired.
754 C JAC need not provide df/dy exactly. A crude
755 C approximation (possibly with greater sparsity) will do.
756 C In any case, PDJ is preset to zero by the solver,
757 C so that only the nonzero elements need be loaded by JAC.
758 C Calls to JAC are made with J = 1,...,NEQ, in that order, and
759 C each such set of calls is preceded by a call to F with the
760 C same arguments NEQ, T, and Y. Thus to gain some efficiency,
761 C intermediate quantities shared by both calculations may be
762 C saved in a user Common block by F and not recomputed by JAC,
763 C if desired. JAC must not alter its input arguments.
764 C JAC must be declared External in the calling program.
765 C Subroutine JAC may access user-defined quantities in
766 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
767 C (dimensioned in JAC) and/or Y has length exceeding NEQ(1).
768 C See the descriptions of NEQ and Y above.
769 C
770 C MF = the method flag. Used only for input.
771 C MF has three decimal digits-- MOSS, METH, MITER--
772 C MF = 100*MOSS + 10*METH + MITER.
773 C MOSS indicates the method to be used to obtain the sparsity
774 C structure of the Jacobian matrix if MITER = 1 or 2:
775 C MOSS = 0 means the user has supplied IA and JA
776 C (see descriptions under IWORK above).
777 C MOSS = 1 means the user has supplied JAC (see below)
778 C and the structure will be obtained from NEQ
779 C initial calls to JAC.
780 C MOSS = 2 means the structure will be obtained from NEQ+1
781 C initial calls to F.
782 C METH indicates the basic linear multistep method:
783 C METH = 1 means the implicit Adams method.
784 C METH = 2 means the method based on Backward
785 C Differentiation Formulas (BDFs).
786 C MITER indicates the corrector iteration method:
787 C MITER = 0 means functional iteration (no Jacobian matrix
788 C is involved).
789 C MITER = 1 means chord iteration with a user-supplied
790 C sparse Jacobian, given by Subroutine JAC.
791 C MITER = 2 means chord iteration with an internally
792 C generated (difference quotient) sparse Jacobian
793 C (using NGP extra calls to F per df/dy value,
794 C where NGP is an optional output described below.)
795 C MITER = 3 means chord iteration with an internally
796 C generated diagonal Jacobian approximation
797 C (using 1 extra call to F per df/dy evaluation).
798 C If MITER = 1 or MOSS = 1, the user must supply a Subroutine
799 C JAC (the name is arbitrary) as described above under JAC.
800 C Otherwise, a dummy argument can be used.
801 C
802 C The standard choices for MF are:
803 C MF = 10 for a nonstiff problem,
804 C MF = 21 or 22 for a stiff problem with IA/JA supplied
805 C (21 if JAC is supplied, 22 if not),
806 C MF = 121 for a stiff problem with JAC supplied,
807 C but not IA/JA,
808 C MF = 222 for a stiff problem with neither IA/JA nor
809 C JAC supplied.
810 C The sparseness structure can be changed during the
811 C problem by making a call to DLSODES with ISTATE = 3.
812 C-----------------------------------------------------------------------
813 C Optional Inputs.
814 C
815 C The following is a list of the optional inputs provided for in the
816 C call sequence. (See also Part 2.) For each such input variable,
817 C this table lists its name as used in this documentation, its
818 C location in the call sequence, its meaning, and the default value.
819 C The use of any of these inputs requires IOPT = 1, and in that
820 C case all of these inputs are examined. A value of zero for any
821 C of these optional inputs will cause the default value to be used.
822 C Thus to use a subset of the optional inputs, simply preload
823 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
824 C then set those of interest to nonzero values.
825 C
826 C Name Location Meaning and Default Value
827 C
828 C H0 RWORK(5) the step size to be attempted on the first step.
829 C The default value is determined by the solver.
830 C
831 C HMAX RWORK(6) the maximum absolute step size allowed.
832 C The default value is infinite.
833 C
834 C HMIN RWORK(7) the minimum absolute step size allowed.
835 C The default value is 0. (This lower bound is not
836 C enforced on the final step before reaching TCRIT
837 C when ITASK = 4 or 5.)
838 C
839 C SETH RWORK(8) the element threshhold for sparsity determination
840 C when MOSS = 1 or 2. If the absolute value of
841 C an estimated Jacobian element is .le. SETH, it
842 C will be assumed to be absent in the structure.
843 C The default value of SETH is 0.
844 C
845 C MAXORD IWORK(5) the maximum order to be allowed. The default
846 C value is 12 if METH = 1, and 5 if METH = 2.
847 C If MAXORD exceeds the default value, it will
848 C be reduced to the default value.
849 C If MAXORD is changed during the problem, it may
850 C cause the current order to be reduced.
851 C
852 C MXSTEP IWORK(6) maximum number of (internally defined) steps
853 C allowed during one call to the solver.
854 C The default value is 500.
855 C
856 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
857 C warning that T + H = T on a step (H = step size).
858 C This must be positive to result in a non-default
859 C value. The default value is 10.
860 C-----------------------------------------------------------------------
861 C Optional Outputs.
862 C
863 C As optional additional output from DLSODES, the variables listed
864 C below are quantities related to the performance of DLSODES
865 C which are available to the user. These are communicated by way of
866 C the work arrays, but also have internal mnemonic names as shown.
867 C Except where stated otherwise, all of these outputs are defined
868 C on any successful return from DLSODES, and on any return with
869 C ISTATE = -1, -2, -4, -5, or -6. On an illegal input return
870 C (ISTATE = -3), they will be unchanged from their existing values
871 C (if any), except possibly for TOLSF, LENRW, and LENIW.
872 C On any error return, outputs relevant to the error will be defined,
873 C as noted below.
874 C
875 C Name Location Meaning
876 C
877 C HU RWORK(11) the step size in t last used (successfully).
878 C
879 C HCUR RWORK(12) the step size to be attempted on the next step.
880 C
881 C TCUR RWORK(13) the current value of the independent variable
882 C which the solver has actually reached, i.e. the
883 C current internal mesh point in t. On output, TCUR
884 C will always be at least as far as the argument
885 C T, but may be farther (if interpolation was done).
886 C
887 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
888 C computed when a request for too much accuracy was
889 C detected (ISTATE = -3 if detected at the start of
890 C the problem, ISTATE = -2 otherwise). If ITOL is
891 C left unaltered but RTOL and ATOL are uniformly
892 C scaled up by a factor of TOLSF for the next call,
893 C then the solver is deemed likely to succeed.
894 C (The user may also ignore TOLSF and alter the
895 C tolerance parameters in any other way appropriate.)
896 C
897 C NST IWORK(11) the number of steps taken for the problem so far.
898 C
899 C NFE IWORK(12) the number of f evaluations for the problem so far,
900 C excluding those for structure determination
901 C (MOSS = 2).
902 C
903 C NJE IWORK(13) the number of Jacobian evaluations for the problem
904 C so far, excluding those for structure determination
905 C (MOSS = 1).
906 C
907 C NQU IWORK(14) the method order last used (successfully).
908 C
909 C NQCUR IWORK(15) the order to be attempted on the next step.
910 C
911 C IMXER IWORK(16) the index of the component of largest magnitude in
912 C the weighted local error vector ( E(i)/EWT(i) ),
913 C on an error return with ISTATE = -4 or -5.
914 C
915 C LENRW IWORK(17) the length of RWORK actually required.
916 C This is defined on normal returns and on an illegal
917 C input return for insufficient storage.
918 C
919 C LENIW IWORK(18) the length of IWORK actually required.
920 C This is defined on normal returns and on an illegal
921 C input return for insufficient storage.
922 C
923 C NNZ IWORK(19) the number of nonzero elements in the Jacobian
924 C matrix, including the diagonal (MITER = 1 or 2).
925 C (This may differ from that given by IA(NEQ+1)-1
926 C if MOSS = 0, because of added diagonal entries.)
927 C
928 C NGP IWORK(20) the number of groups of column indices, used in
929 C difference quotient Jacobian aproximations if
930 C MITER = 2. This is also the number of extra f
931 C evaluations needed for each Jacobian evaluation.
932 C
933 C NLU IWORK(21) the number of sparse LU decompositions for the
934 C problem so far.
935 C
936 C LYH IWORK(22) the base address in RWORK of the history array YH,
937 C described below in this list.
938 C
939 C IPIAN IWORK(23) the base address of the structure descriptor array
940 C IAN, described below in this list.
941 C
942 C IPJAN IWORK(24) the base address of the structure descriptor array
943 C JAN, described below in this list.
944 C
945 C NZL IWORK(25) the number of nonzero elements in the strict lower
946 C triangle of the LU factorization used in the chord
947 C iteration (MITER = 1 or 2).
948 C
949 C NZU IWORK(26) the number of nonzero elements in the strict upper
950 C triangle of the LU factorization used in the chord
951 C iteration (MITER = 1 or 2).
952 C The total number of nonzeros in the factorization
953 C is therefore NZL + NZU + NEQ.
954 C
955 C The following four arrays are segments of the RWORK array which
956 C may also be of interest to the user as optional outputs.
957 C For each array, the table below gives its internal name,
958 C its base address, and its description.
959 C For YH and ACOR, the base addresses are in RWORK (a real array).
960 C The integer arrays IAN and JAN are to be obtained by declaring an
961 C integer array IWK and identifying IWK(1) with RWORK(21), using either
962 C an equivalence statement or a subroutine call. Then the base
963 C addresses IPIAN (of IAN) and IPJAN (of JAN) in IWK are to be obtained
964 C as optional outputs IWORK(23) and IWORK(24), respectively.
965 C Thus IAN(1) is IWK(IPIAN), etc.
966 C
967 C Name Base Address Description
968 C
969 C IAN IPIAN (in IWK) structure descriptor array of size NEQ + 1.
970 C JAN IPJAN (in IWK) structure descriptor array of size NNZ.
971 C (see above) IAN and JAN together describe the sparsity
972 C structure of the Jacobian matrix, as used by
973 C DLSODES when MITER = 1 or 2.
974 C JAN contains the row indices of the nonzero
975 C locations, reading in columnwise order, and
976 C IAN contains the starting locations in JAN of
977 C the descriptions of columns 1,...,NEQ, in
978 C that order, with IAN(1) = 1. Thus for each
979 C j = 1,...,NEQ, the row indices i of the
980 C nonzero locations in column j are
981 C i = JAN(k), IAN(j) .le. k .lt. IAN(j+1).
982 C Note that IAN(NEQ+1) = NNZ + 1.
983 C (If MOSS = 0, IAN/JAN may differ from the
984 C input IA/JA because of a different ordering
985 C in each column, and added diagonal entries.)
986 C
987 C YH LYH the Nordsieck history array, of size NYH by
988 C (optional (NQCUR + 1), where NYH is the initial value
989 C output) of NEQ. For j = 0,1,...,NQCUR, column j+1
990 C of YH contains HCUR**j/factorial(j) times
991 C the j-th derivative of the interpolating
992 C polynomial currently representing the solution,
993 C evaluated at t = TCUR. The base address LYH
994 C is another optional output, listed above.
995 C
996 C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
997 C corrections on each step, scaled on output
998 C to represent the estimated local error in y
999 C on the last step. This is the vector E in
1000 C the description of the error control. It is
1001 C defined only on a successful return from
1002 C DLSODES.
1003 C
1004 C-----------------------------------------------------------------------
1005 C Part 2. Other Routines Callable.
1006 C
1007 C The following are optional calls which the user may make to
1008 C gain additional capabilities in conjunction with DLSODES.
1009 C (The routines XSETUN and XSETF are designed to conform to the
1010 C SLATEC error handling package.)
1011 C
1012 C Form of Call Function
1013 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
1014 C output of messages from DLSODES, if
1015 C the default is not desired.
1016 C The default value of LUN is 6.
1017 C
1018 C CALL XSETF(MFLAG) Set a flag to control the printing of
1019 C messages by DLSODES.
1020 C MFLAG = 0 means do not print. (Danger:
1021 C This risks losing valuable information.)
1022 C MFLAG = 1 means print (the default).
1023 C
1024 C Either of the above calls may be made at
1025 C any time and will take effect immediately.
1026 C
1027 C CALL DSRCMS(RSAV,ISAV,JOB) saves and restores the contents of
1028 C the internal Common blocks used by
1029 C DLSODES (see Part 3 below).
1030 C RSAV must be a real array of length 224
1031 C or more, and ISAV must be an integer
1032 C array of length 71 or more.
1033 C JOB=1 means save Common into RSAV/ISAV.
1034 C JOB=2 means restore Common from RSAV/ISAV.
1035 C DSRCMS is useful if one is
1036 C interrupting a run and restarting
1037 C later, or alternating between two or
1038 C more problems solved with DLSODES.
1039 C
1040 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
1041 C (see below) orders, at a specified point t, if
1042 C desired. It may be called only after
1043 C a successful return from DLSODES.
1044 C
1045 C The detailed instructions for using DINTDY are as follows.
1046 C The form of the call is:
1047 C
1048 C LYH = IWORK(22)
1049 C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
1050 C
1051 C The input parameters are:
1052 C
1053 C T = value of independent variable where answers are desired
1054 C (normally the same as the T last returned by DLSODES).
1055 C For valid results, T must lie between TCUR - HU and TCUR.
1056 C (See optional outputs for TCUR and HU.)
1057 C K = integer order of the derivative desired. K must satisfy
1058 C 0 .le. K .le. NQCUR, where NQCUR is the current order
1059 C (See optional outputs). The capability corresponding
1060 C to K = 0, i.e. computing y(T), is already provided
1061 C by DLSODES directly. Since NQCUR .ge. 1, the first
1062 C derivative dy/dt is always available with DINTDY.
1063 C LYH = the base address of the history array YH, obtained
1064 C as an optional output as shown above.
1065 C NYH = column length of YH, equal to the initial value of NEQ.
1066 C
1067 C The output parameters are:
1068 C
1069 C DKY = a real array of length NEQ containing the computed value
1070 C of the K-th derivative of y(t).
1071 C IFLAG = integer flag, returned as 0 if K and T were legal,
1072 C -1 if K was illegal, and -2 if T was illegal.
1073 C On an error return, a message is also written.
1074 C-----------------------------------------------------------------------
1075 C Part 3. Common Blocks.
1076 C
1077 C If DLSODES is to be used in an overlay situation, the user
1078 C must declare, in the primary overlay, the variables in:
1079 C (1) the call sequence to DLSODES, and
1080 C (2) the two internal Common blocks
1081 C /DLS001/ of length 255 (218 double precision words
1082 C followed by 37 integer words),
1083 C /DLSS01/ of length 40 (6 double precision words
1084 C followed by 34 integer words),
1085 C
1086 C If DLSODES is used on a system in which the contents of internal
1087 C Common blocks are not preserved between calls, the user should
1088 C declare the above Common blocks in the calling program to insure
1089 C that their contents are preserved.
1090 C
1091 C If the solution of a given problem by DLSODES is to be interrupted
1092 C and then later continued, such as when restarting an interrupted run
1093 C or alternating between two or more problems, the user should save,
1094 C following the return from the last DLSODES call prior to the
1095 C interruption, the contents of the call sequence variables and the
1096 C internal Common blocks, and later restore these values before the
1097 C next DLSODES call for that problem. To save and restore the Common
1098 C blocks, use Subroutine DSRCMS (see Part 2 above).
1099 C
1100 C-----------------------------------------------------------------------
1101 C Part 4. Optionally Replaceable Solver Routines.
1102 C
1103 C Below are descriptions of two routines in the DLSODES package which
1104 C relate to the measurement of errors. Either routine can be
1105 C replaced by a user-supplied version, if desired. However, since such
1106 C a replacement may have a major impact on performance, it should be
1107 C done only when absolutely necessary, and only with great caution.
1108 C (Note: The means by which the package version of a routine is
1109 C superseded by the user's version may be system-dependent.)
1110 C
1111 C (a) DEWSET.
1112 C The following subroutine is called just before each internal
1113 C integration step, and sets the array of error weights, EWT, as
1114 C described under ITOL/RTOL/ATOL above:
1115 C Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1116 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODES call sequence,
1117 C YCUR contains the current dependent variable vector, and
1118 C EWT is the array of weights set by DEWSET.
1119 C
1120 C If the user supplies this subroutine, it must return in EWT(i)
1121 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1122 C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
1123 C routine (see below), and also used by DLSODES in the computation
1124 C of the optional output IMXER, the diagonal Jacobian approximation,
1125 C and the increments for difference quotient Jacobians.
1126 C
1127 C In the user-supplied version of DEWSET, it may be desirable to use
1128 C the current values of derivatives of y. Derivatives up to order NQ
1129 C are available from the history array YH, described above under
1130 C optional outputs. In DEWSET, YH is identical to the YCUR array,
1131 C extended to NQ + 1 columns with a column length of NYH and scale
1132 C factors of H**j/factorial(j). On the first call for the problem,
1133 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1134 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
1135 C can be obtained by including in DEWSET the statements:
1136 C DOUBLE PRECISION RLS
1137 C COMMON /DLS001/ RLS(218),ILS(37)
1138 C NQ = ILS(33)
1139 C NST = ILS(34)
1140 C H = RLS(212)
1141 C Thus, for example, the current value of dy/dt can be obtained as
1142 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
1143 C unnecessary when NST = 0).
1144 C
1145 C (b) DVNORM.
1146 C The following is a real function routine which computes the weighted
1147 C root-mean-square norm of a vector v:
1148 C D = DVNORM (N, V, W)
1149 C where
1150 C N = the length of the vector,
1151 C V = real array of length N containing the vector,
1152 C W = real array of length N containing weights,
1153 C D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).
1154 C DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
1155 C EWT is as set by Subroutine DEWSET.
1156 C
1157 C If the user supplies this function, it should return a non-negative
1158 C value of DVNORM suitable for use in the error control in DLSODES.
1159 C None of the arguments should be altered by DVNORM.
1160 C For example, a user-supplied DVNORM routine might:
1161 C -substitute a max-norm of (V(i)*W(i)) for the RMS-norm, or
1162 C -ignore some components of V in the norm, with the effect of
1163 C suppressing the error control on those components of y.
1164 C-----------------------------------------------------------------------
1165 C
1166 C***REVISION HISTORY (YYYYMMDD)
1167 C 19810120 DATE WRITTEN
1168 C 19820315 Upgraded MDI in ODRV package: operates on M + M-transpose.
1169 C 19820426 Numerous revisions in use of work arrays;
1170 C use wordlength ratio LENRAT; added IPISP & LRAT to Common;
1171 C added optional outputs IPIAN/IPJAN;
1172 C numerous corrections to comments.
1173 C 19830503 Added routine CNTNZU; added NZL and NZU to /LSS001/;
1174 C changed ADJLR call logic; added optional outputs NZL & NZU;
1175 C revised counter initializations; revised PREP stmt. numbers;
1176 C corrections to comments throughout.
1177 C 19870320 Corrected jump on test of umax in CDRV routine;
1178 C added ISTATE = -7 return.
1179 C 19870330 Major update: corrected comments throughout;
1180 C removed TRET from Common; rewrote EWSET with 4 loops;
1181 C fixed t test in INTDY; added Cray directives in STODE;
1182 C in STODE, fixed DELP init. and logic around PJAC call;
1183 C combined routines to save/restore Common;
1184 C passed LEVEL = 0 in error message calls (except run abort).
1185 C 20010425 Major update: convert source lines to upper case;
1186 C added *DECK lines; changed from 1 to * in dummy dimensions;
1187 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
1188 C renamed routines for uniqueness across single/double prec.;
1189 C converted intrinsic names to generic form;
1190 C removed ILLIN and NTREP (data loaded) from Common;
1191 C removed all 'own' variables from Common;
1192 C changed error messages to quoted strings;
1193 C replaced XERRWV/XERRWD with 1993 revised version;
1194 C converted prologues, comments, error messages to mixed case;
1195 C converted arithmetic IF statements to logical IF statements;
1196 C numerous corrections to prologues and internal comments.
1197 C 20010507 Converted single precision source to double precision.
1198 C 20020502 Corrected declarations in descriptions of user routines.
1199 C 20031105 Restored 'own' variables to Common blocks, to enable
1200 C interrupt/restart feature.
1201 C 20031112 Added SAVE statements for data-loaded constants.
1202 C
1203 C-----------------------------------------------------------------------
1204 C Other routines in the DLSODES package.
1205 C
1206 C In addition to Subroutine DLSODES, the DLSODES package includes the
1207 C following subroutines and function routines:
1208 C DIPREP acts as an iterface between DLSODES and DPREP, and also does
1209 C adjusting of work space pointers and work arrays.
1210 C DPREP is called by DIPREP to compute sparsity and do sparse matrix
1211 C preprocessing if MITER = 1 or 2.
1212 C JGROUP is called by DPREP to compute groups of Jacobian column
1213 C indices for use when MITER = 2.
1214 C ADJLR adjusts the length of required sparse matrix work space.
1215 C It is called by DPREP.
1216 C CNTNZU is called by DPREP and counts the nonzero elements in the
1217 C strict upper triangle of J + J-transpose, where J = df/dy.
1218 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1219 C DSTODE is the core integrator, which does one step of the
1220 C integration and the associated error control.
1221 C DCFODE sets all method coefficients and test constants.
1222 C DPRJS computes and preprocesses the Jacobian matrix J = df/dy
1223 C and the Newton iteration matrix P = I - h*l0*J.
1224 C DSOLSS manages solution of linear system in chord iteration.
1225 C DEWSET sets the error weight vector EWT before each step.
1226 C DVNORM computes the weighted RMS-norm of a vector.
1227 C DSRCMS is a user-callable routine to save and restore
1228 C the contents of the internal Common blocks.
1229 C ODRV constructs a reordering of the rows and columns of
1230 C a matrix by the minimum degree algorithm. ODRV is a
1231 C driver routine which calls Subroutines MD, MDI, MDM,
1232 C MDP, MDU, and SRO. See Ref. 2 for details. (The ODRV
1233 C module has been modified since Ref. 2, however.)
1234 C CDRV performs reordering, symbolic factorization, numerical
1235 C factorization, or linear system solution operations,
1236 C depending on a path argument ipath. CDRV is a
1237 C driver routine which calls Subroutines NROC, NSFC,
1238 C NNFC, NNSC, and NNTC. See Ref. 3 for details.
1239 C DLSODES uses CDRV to solve linear systems in which the
1240 C coefficient matrix is P = I - con*J, where I is the
1241 C identity, con is a scalar, and J is an approximation to
1242 C the Jacobian df/dy. Because CDRV deals with rowwise
1243 C sparsity descriptions, CDRV works with P-transpose, not P.
1244 C DUMACH computes the unit roundoff in a machine-independent manner.
1245 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
1246 C error messages and warnings. XERRWD is machine-dependent.
1247 C Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1248 C All the others are subroutines.
1249 C
1250 C-----------------------------------------------------------------------
1270 LOGICAL IHIT
1271 CHARACTER*60 MSG
1273 C-----------------------------------------------------------------------
1274 C The following two internal Common blocks contain
1275 C (a) variables which are local to any subroutine but whose values must
1276 C be preserved between calls to the routine ("own" variables), and
1277 C (b) variables which are communicated between subroutines.
1278 C The block DLS001 is declared in subroutines DLSODES, DIPREP, DPREP,
1279 C DINTDY, DSTODE, DPRJS, and DSOLSS.
1280 C The block DLSS01 is declared in subroutines DLSODES, DIPREP, DPREP,
1281 C DPRJS, and DSOLSS.
1282 C Groups of variables are replaced by dummy arrays in the Common
1283 C declarations in routines where those variables are not used.
1284 C-----------------------------------------------------------------------
1291 C
1297 C
1299 C-----------------------------------------------------------------------
1300 C In the Data statement below, set LENRAT equal to the ratio of
1301 C the wordlength for a real number to that for an integer. Usually,
1302 C LENRAT = 1 for single precision and 2 for double precision. If the
1303 C true ratio is not an integer, use the next smaller integer (.ge. 1).
1304 C-----------------------------------------------------------------------
1306 C-----------------------------------------------------------------------
1307 C Block A.
1308 C This code block is executed on every call.
1309 C It tests ISTATE and ITASK for legality and branches appropriately.
1310 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1311 C not yet been done, an error return occurs.
1312 C If ISTATE = 1 and TOUT = T, return immediately.
1313 C-----------------------------------------------------------------------
1322 C-----------------------------------------------------------------------
1323 C Block B.
1324 C The next code block is executed for the initial call (ISTATE = 1),
1325 C or for a continuation call with parameter changes (ISTATE = 3).
1326 C It contains checking of all inputs and various initializations.
1327 C If ISTATE = 1, the final setting of work space pointers, the matrix
1328 C preprocessing, and other initializations are done in Block C.
1329 C
1330 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1331 C MF, ML, and MU.
1332 C-----------------------------------------------------------------------
1347 C Next process and check the optional inputs. --------------------------
1350 MXSTEP = MXSTP0
1351 MXHNIL = MXHNL0
1378 C Check RTOL and ATOL for legality. ------------------------------------
1387 C-----------------------------------------------------------------------
1388 C Compute required work array lengths, as far as possible, and test
1389 C these against LRW and LIW. Then set tentative pointers for work
1390 C arrays. Pointers to RWORK/IWORK segments are named by prefixing L to
1391 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1392 C Segments of RWORK (in order) are denoted WM, YH, SAVF, EWT, ACOR.
1393 C If MITER = 1 or 2, the required length of the matrix work space WM
1394 C is not yet known, and so a crude minimum value is used for the
1395 C initial tests of LRW and LIW, and YH is temporarily stored as far
1396 C to the right in RWORK as possible, to leave the maximum amount
1397 C of space for WM for matrix preprocessing. Thus if MITER = 1 or 2
1398 C and MOSS .ne. 2, some of the segments of RWORK are temporarily
1399 C omitted, as they are not needed in the preprocessing. These
1400 C omitted segments are: ACOR if ISTATE = 1, EWT and ACOR if ISTATE = 3
1401 C and MOSS = 1, and SAVF, EWT, and ACOR if ISTATE = 3 and MOSS = 0.
1402 C-----------------------------------------------------------------------
1403 LRAT = LENRAT
1430 LENYHM = NCOLM*NYH
1431 LENYHT = LENYH
1437 LWTEM = LWMIN
1439 LENWK = LWTEM
1444 ISTATC = ISTATE
1446 C-----------------------------------------------------------------------
1447 C ISTATE = 3. Move YH to its new location.
1448 C Note that only the part of YH needed for the next step, namely
1449 C MIN(NQ+1,MAXORD+2) columns, is actually moved.
1450 C A temporary error weight array EWT is loaded if MOSS = 2.
1451 C Sparse matrix processing is done in DIPREP/DPREP if MITER = 1 or 2.
1452 C If MAXORD was reduced below NQ, then the pointers are finally set
1453 C so that SAVF is identical to YH(*,MAXORD+2).
1454 C-----------------------------------------------------------------------
1457 C Move YH. Move right if LYHD < 0; move left if LYHD > 0. -------------
1462 ENDIF
1466 ENDIF
1471 C Temporarily load EWT if MITER = 1 or 2 and MOSS = 2. -----------------
1477 C DIPREP and DPREP do sparse matrix preprocessing if MITER = 1 or 2. ---
1490 C Set flag to signal parameter changes to DSTODE. ----------------------
1493 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1500 C-----------------------------------------------------------------------
1501 C Block C.
1502 C The next block is for the initial call only (ISTATE = 1).
1503 C It contains all remaining initializations, the initial call to F,
1504 C the sparse matrix preprocessing (MITER = 1 or 2), and the
1505 C calculation of the initial step size.
1506 C The error weights in EWT are inverted after being loaded.
1507 C-----------------------------------------------------------------------
1509 LYH = LYHN
1511 TN = T
1518 C Load the initial value vector in YH. ---------------------------------
1521 C Initial call to F. (LF0 points to YH(*,2).) -------------------------
1525 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1531 C DIPREP and DPREP do sparse matrix preprocessing if MITER = 1 or 2. ---
1542 C Check TCRIT for legality (ITASK = 4 or 5). ---------------------------
1549 C Initialize all remaining parameters. ---------------------------------
1568 C-----------------------------------------------------------------------
1569 C The coding below computes the step size, H0, to be attempted on the
1570 C first step, unless the user has supplied a value for this.
1571 C First check that TOUT - T differs significantly from zero.
1572 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1573 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1574 C so as to be between 100*UROUND and 1.0E-3.
1575 C Then the computed value H0 is given by..
1576 C NEQ
1577 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * Sum ( f(i)/ywt(i) )**2 )
1578 C 1
1579 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1580 C f(i) = i-th component of initial value of f,
1581 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1582 C The sign of H0 is inferred from the initial values of TOUT and T.
1583 C ABS(H0) is made .le. ABS(TOUT-T) in any case.
1584 C-----------------------------------------------------------------------
1608 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1611 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1612 H = H0
1616 C-----------------------------------------------------------------------
1617 C Block D.
1618 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1619 C and is to check stop conditions before taking a step.
1620 C-----------------------------------------------------------------------
1626 T = TOUT
1638 T = TOUT
1649 C-----------------------------------------------------------------------
1650 C Block E.
1651 C The next block is normally executed for all calls and contains
1652 C the call to the one-step core integrator DSTODE.
1653 C
1654 C This is a looping point for the integration steps.
1655 C
1656 C First check for too many steps being taken, update EWT (if not at
1657 C start of problem), check for too much accuracy being requested, and
1658 C check for H below the roundoff level in T.
1659 C-----------------------------------------------------------------------
1686 C-----------------------------------------------------------------------
1687 C CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,WM,F,JAC,DPRJS,DSOLSS)
1688 C-----------------------------------------------------------------------
1694 C-----------------------------------------------------------------------
1695 C Block F.
1696 C The following block handles the case of a successful return from the
1697 C core integrator (KFLAG = 0). Test for stop conditions.
1698 C-----------------------------------------------------------------------
1701 C ITASK = 1. if TOUT has been reached, interpolate. -------------------
1704 T = TOUT
1706 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1709 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1712 T = TOUT
1722 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1725 C-----------------------------------------------------------------------
1726 C Block G.
1727 C The following block handles all successful returns from DLSODES.
1728 C If ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1729 C ISTATE is set to 2, and the optional outputs are loaded into the
1730 C work arrays before returning.
1731 C-----------------------------------------------------------------------
1734 T = TN
1751 RETURN
1752 C-----------------------------------------------------------------------
1753 C Block H.
1754 C The following block handles all unsuccessful returns other than
1755 C those for illegal input. First the error message routine is called.
1756 C If there was an error test or convergence test failure, IMXER is set.
1757 C Then Y is loaded from YH and T is set to TN.
1758 C The optional outputs are loaded into the work arrays before returning.
1759 C-----------------------------------------------------------------------
1760 C The maximum number of steps was taken before reaching TOUT. ----------
1767 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1773 C Too much accuracy requested for machine precision. -------------------
1781 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1788 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1797 C KFLAG = -3. Fatal error flag returned by DPRJS or DSOLSS (CDRV). ----
1806 C Compute IMXER if relevant. -------------------------------------------
1812 BIG = SIZE
1813 IMXER = I
1816 C Set Y vector, T, and optional outputs. -------------------------------
1819 T = TN
1833 RETURN
1834 C-----------------------------------------------------------------------
1835 C Block I.
1836 C The following block handles all error returns due to illegal input
1837 C (ISTATE = -3), as detected before calling the core integrator.
1838 C First the error message routine is called. If the illegal input
1839 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1840 C-----------------------------------------------------------------------
1966 ENDIF
1970 ENDIF
1971 C
1973 RETURN
1974 C
1977 RETURN
1978 C----------------------- End of Subroutine DLSODES ---------------------
1979 END