| blob | f2926a7201ed0d148a3b2c24466684fc2d34ba9e |
1 *DECK DLSODIS
9 C-----------------------------------------------------------------------
10 C This is the 18 November 2003 version of
11 C DLSODIS: Livermore Solver for Ordinary Differential equations
12 C (Implicit form) with general Sparse Jacobian matrices.
13 C
14 C This version is in double precision.
15 C
16 C DLSODIS solves the initial value problem for linearly implicit
17 C systems of first order ODEs,
18 C A(t,y) * dy/dt = g(t,y) , where A(t,y) is a square matrix,
19 C or, in component form,
20 C ( a * ( dy / dt )) + ... + ( a * ( dy / dt )) =
21 C i,1 1 i,NEQ NEQ
22 C
23 C = g ( t, y , y ,..., y ) ( i = 1,...,NEQ )
24 C i 1 2 NEQ
25 C
26 C If A is singular, this is a differential-algebraic system.
27 C
28 C DLSODIS is a variant version of the DLSODI package, and is intended
29 C for stiff problems in which the matrix A and the Jacobian matrix
30 C d(g - A*s)/dy have arbitrary sparse structures.
31 C
32 C Authors: Alan C. Hindmarsh
33 C Center for Applied Scientific Computing, L-561
34 C Lawrence Livermore National Laboratory
35 C Livermore, CA 94551
36 C and
37 C Sheila Balsdon
38 C Zycor, Inc.
39 C Austin, TX 78741
40 C-----------------------------------------------------------------------
41 C References:
42 C 1. M. K. Seager and S. Balsdon, LSODIS, A Sparse Implicit
43 C ODE Solver, in Proceedings of the IMACS 10th World Congress,
44 C Montreal, August 8-13, 1982.
45 C
46 C 2. Alan C. Hindmarsh, LSODE and LSODI, Two New Initial Value
47 C Ordinary Differential Equation Solvers,
48 C ACM-SIGNUM Newsletter, vol. 15, no. 4 (1980), pp. 10-11.
49 C
50 C 3. S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman,
51 C Yale Sparse Matrix Package: I. The Symmetric Codes,
52 C Int. J. Num. Meth. Eng., vol. 18 (1982), pp. 1145-1151.
53 C
54 C 4. S. C. Eisenstat, M. C. Gursky, M. H. Schultz, and A. H. Sherman,
55 C Yale Sparse Matrix Package: II. The Nonsymmetric Codes,
56 C Research Report No. 114, Dept. of Computer Sciences, Yale
57 C University, 1977.
58 C-----------------------------------------------------------------------
59 C Summary of Usage.
60 C
61 C Communication between the user and the DLSODIS package, for normal
62 C situations, is summarized here. This summary describes only a subset
63 C of the full set of options available. See the full description for
64 C details, including optional communication, nonstandard options,
65 C and instructions for special situations. See also the example
66 C problem (with program and output) following this summary.
67 C
68 C A. First, provide a subroutine of the form:
69 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
70 C DOUBLE PRECISION T, Y(*), S(*), R(*)
71 C which computes the residual function
72 C r = g(t,y) - A(t,y) * s ,
73 C as a function of t and the vectors y and s. (s is an internally
74 C generated approximation to dy/dt.) The arrays Y and S are inputs
75 C to the RES routine and should not be altered. The residual
76 C vector is to be stored in the array R. The argument IRES should be
77 C ignored for casual use of DLSODIS. (For uses of IRES, see the
78 C paragraph on RES in the full description below.)
79 C
80 C B. DLSODIS must deal internally with the matrices A and dr/dy, where
81 C r is the residual function defined above. DLSODIS generates a linear
82 C combination of these two matrices in sparse form.
83 C The matrix structure is communicated by a method flag, MF:
84 C MF = 21 or 22 when the user provides the structures of
85 C matrix A and dr/dy,
86 C MF = 121 or 222 when the user does not provide structure
87 C information, and
88 C MF = 321 or 422 when the user provides the structure
89 C of matrix A.
90 C
91 C C. You must also provide a subroutine of the form:
92 C SUBROUTINE ADDA (NEQ, T, Y, J, IAN, JAN, P)
93 C DOUBLE PRECISION T, Y(*), P(*)
94 C INTEGER IAN(*), JAN(*)
95 C which adds the matrix A = A(t,y) to the contents of the array P.
96 C NEQ, T, Y, and J are input arguments and should not be altered.
97 C This routine should add the J-th column of matrix A to the array
98 C P (of length NEQ). I.e. add A(i,J) to P(i) for all relevant
99 C values of i. The arguments IAN and JAN should be ignored for normal
100 C situations. DLSODIS will call the ADDA routine with J = 1,2,...,NEQ.
101 C
102 C D. For the sake of efficiency, you are encouraged to supply the
103 C Jacobian matrix dr/dy in closed form, where r = g(t,y) - A(t,y)*s
104 C (s = a fixed vector) as above. If dr/dy is being supplied,
105 C use MF = 21, 121, or 321, and provide a subroutine of the form:
106 C SUBROUTINE JAC (NEQ, T, Y, S, J, IAN, JAN, PDJ)
107 C DOUBLE PRECISION T, Y(*), S(*), PDJ(*)
108 C INTEGER IAN(*), JAN(*)
109 C which computes dr/dy as a function of t, y, and s. Here NEQ, T, Y, S,
110 C and J are input arguments, and the JAC routine is to load the array
111 C PDJ (of length NEQ) with the J-th column of dr/dy. I.e. load PDJ(i)
112 C with dr(i)/dy(J) for all relevant values of i. The arguments IAN and
113 C JAN should be ignored for normal situations. DLSODIS will call the
114 C JAC routine with J = 1,2,...,NEQ.
115 C Only nonzero elements need be loaded. A crude approximation
116 C to dr/dy, possibly with fewer nonzero elememts, will suffice.
117 C Note that if A is independent of y (or this dependence
118 C is weak enough to be ignored) then JAC is to compute dg/dy.
119 C If it is not feasible to provide a JAC routine, use
120 C MF = 22, 222, or 422 and DLSODIS will compute an approximate
121 C Jacobian internally by difference quotients.
122 C
123 C E. Next decide whether or not to provide the initial value of the
124 C derivative vector dy/dt. If the initial value of A(t,y) is
125 C nonsingular (and not too ill-conditioned), you may let DLSODIS compute
126 C this vector (ISTATE = 0). (DLSODIS will solve the system A*s = g for
127 C s, with initial values of A and g.) If A(t,y) is initially
128 C singular, then the system is a differential-algebraic system, and
129 C you must make use of the particular form of the system to compute the
130 C initial values of y and dy/dt. In that case, use ISTATE = 1 and
131 C load the initial value of dy/dt into the array YDOTI.
132 C The input array YDOTI and the initial Y array must be consistent with
133 C the equations A*dy/dt = g. This implies that the initial residual
134 C r = g(t,y) - A(t,y)*YDOTI must be approximately zero.
135 C
136 C F. Write a main program which calls Subroutine DLSODIS once for
137 C each point at which answers are desired. This should also provide
138 C for possible use of logical unit 6 for output of error messages by
139 C DLSODIS. On the first call to DLSODIS, supply arguments as follows:
140 C RES = name of user subroutine for residual function r.
141 C ADDA = name of user subroutine for computing and adding A(t,y).
142 C JAC = name of user subroutine for Jacobian matrix dr/dy
143 C (MF = 121). If not used, pass a dummy name.
144 C Note: The names for the RES and ADDA routines and (if used) the
145 C JAC routine must be declared External in the calling program.
146 C NEQ = number of scalar equations in the system.
147 C Y = array of initial values, of length NEQ.
148 C YDOTI = array of length NEQ (containing initial dy/dt if ISTATE = 1).
149 C T = the initial value of the independent variable.
150 C TOUT = first point where output is desired (.ne. T).
151 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
152 C RTOL = relative tolerance parameter (scalar).
153 C ATOL = absolute tolerance parameter (scalar or array).
154 C The estimated local error in y(i) will be controlled so as
155 C to be roughly less (in magnitude) than
156 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
157 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
158 C Thus the local error test passes if, in each component,
159 C either the absolute error is less than ATOL (or ATOL(i)),
160 C or the relative error is less than RTOL.
161 C Use RTOL = 0.0 for pure absolute error control, and
162 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
163 C control. Caution: Actual (global) errors may exceed these
164 C local tolerances, so choose them conservatively.
165 C ITASK = 1 for normal computation of output values of y at t = TOUT.
166 C ISTATE = integer flag (input and output). Set ISTATE = 1 if the
167 C initial dy/dt is supplied, and 0 otherwise.
168 C IOPT = 0 to indicate no optional inputs used.
169 C RWORK = real work array of length at least:
170 C 20 + (2 + 1./LENRAT)*NNZ + (11 + 9./LENRAT)*NEQ
171 C where:
172 C NNZ = the number of nonzero elements in the sparse
173 C iteration matrix P = A - con*dr/dy (con = scalar)
174 C (If NNZ is unknown, use an estimate of it.)
175 C LENRAT = the real to integer wordlength ratio (usually 1 in
176 C single precision and 2 in double precision).
177 C In any case, the required size of RWORK cannot generally
178 C be predicted in advance for any value of MF, and the
179 C value above is a rough estimate of a crude lower bound.
180 C Some experimentation with this size may be necessary.
181 C (When known, the correct required length is an optional
182 C output, available in IWORK(17).)
183 C LRW = declared length of RWORK (in user's dimension).
184 C IWORK = integer work array of length at least 30.
185 C LIW = declared length of IWORK (in user's dimension).
186 C MF = method flag. Standard values are:
187 C 121 for a user-supplied sparse Jacobian.
188 C 222 for an internally generated sparse Jacobian.
189 C For other choices of MF, see the paragraph on MF in
190 C the full description below.
191 C Note that the main program must declare arrays Y, YDOTI, RWORK, IWORK,
192 C and possibly ATOL.
193 C
194 C G. The output from the first call, or any call, is:
195 C Y = array of computed values of y(t) vector.
196 C T = corresponding value of independent variable (normally TOUT).
197 C ISTATE = 2 if DLSODIS was successful, negative otherwise.
198 C -1 means excess work done on this call (check all inputs).
199 C -2 means excess accuracy requested (tolerances too small).
200 C -3 means illegal input detected (see printed message).
201 C -4 means repeated error test failures (check all inputs).
202 C -5 means repeated convergence failures (perhaps bad Jacobian
203 C supplied or wrong choice of tolerances).
204 C -6 means error weight became zero during problem. (Solution
205 C component i vanished, and ATOL or ATOL(i) = 0.)
206 C -7 cannot occur in casual use.
207 C -8 means DLSODIS was unable to compute the initial dy/dt.
208 C in casual use, this means A(t,y) is initially singular.
209 C Supply YDOTI and use ISTATE = 1 on the first call.
210 C -9 means a fatal error return flag came from sparse solver
211 C CDRV by way of DPRJIS or DSOLSS. Should never happen.
212 C
213 C A return with ISTATE = -1, -4, or -5, may result from using
214 C an inappropriate sparsity structure, one that is quite
215 C different from the initial structure. Consider calling
216 C DLSODIS again with ISTATE = 3 to force the structure to be
217 C reevaluated. See the full description of ISTATE below.
218 C
219 C If DLSODIS returns ISTATE = -1, -4 or -5, then the output of
220 C DLSODIS also includes YDOTI = array containing residual vector
221 C r = g - A * dy/dt evaluated at the current t, y, and dy/dt.
222 C
223 C H. To continue the integration after a successful return, simply
224 C reset TOUT and call DLSODIS again. No other parameters need be reset.
225 C
226 C-----------------------------------------------------------------------
227 C Example Problem.
228 C
229 C The following is an example problem, with the coding needed
230 C for its solution by DLSODIS. The problem comes from the partial
231 C differential equation (the Burgers equation)
232 C du/dt = - u * du/dx + eta * d**2 u/dx**2, eta = .05,
233 C on -1 .le. x .le. 1. The boundary conditions are periodic:
234 C u(-1,t) = u(1,t) and du/dx(-1,t) = du/dx(1,t)
235 C The initial profile is a square wave,
236 C u = 1 in ABS(x) .lt. .5, u = .5 at ABS(x) = .5, u = 0 elsewhere.
237 C The PDE is discretized in x by a simplified Galerkin method,
238 C using piecewise linear basis functions, on a grid of 40 intervals.
239 C The result is a system A * dy/dt = g(y), of size NEQ = 40,
240 C where y(i) is the approximation to u at x = x(i), with
241 C x(i) = -1 + (i-1)*delx, delx = 2/NEQ = .05.
242 C The individual equations in the system are (in order):
243 C (1/6)dy(NEQ)/dt+(4/6)dy(1)/dt+(1/6)dy(2)/dt
244 C = r4d*(y(NEQ)**2-y(2)**2)+eodsq*(y(2)-2*y(1)+y(NEQ))
245 C for i = 2,3,...,nm1,
246 C (1/6)dy(i-1)/dt+(4/6)dy(i)/dt+(1/6)dy(i+1)/dt
247 C = r4d*(y(i-1)**2-y(i+1)**2)+eodsq*(y(i+1)-2*y(i)+y(i-1))
248 C and finally
249 C (1/6)dy(nm1)/dt+(4/6)dy(NEQ)/dt+(1/6)dy(1)/dt
250 C = r4d*(y(nm1)**2-y(1)**2)+eodsq*(y(1)-2*y(NEQ)+y(nm1))
251 C where r4d = 1/(4*delx), eodsq = eta/delx**2 and nm1 = NEQ-1.
252 C The following coding solves the problem with MF = 121, with output
253 C of solution statistics at t = .1, .2, .3, and .4, and of the
254 C solution vector at t = .4. Optional outputs (run statistics) are
255 C also printed.
256 C
257 C EXTERNAL RESID, ADDASP, JACSP
258 C DOUBLE PRECISION ATOL, RTOL, RW, T, TOUT, Y, YDOTI, R4D, EODSQ, DELX
259 C DIMENSION Y(40), YDOTI(40), RW(1409), IW(30)
260 C COMMON /TEST1/ R4D, EODSQ, NM1
261 C DATA ITOL/1/, RTOL/1.0D-3/, ATOL/1.0D-3/, ITASK/1/, IOPT/0/
262 C DATA NEQ/40/, LRW/1409/, LIW/30/, MF/121/
263 C
264 C DELX = 2.0/NEQ
265 C R4D = 0.25/DELX
266 C EODSQ = 0.05/DELX**2
267 C NM1 = NEQ - 1
268 C DO 10 I = 1,NEQ
269 C 10 Y(I) = 0.0
270 C Y(11) = 0.5
271 C DO 15 I = 12,30
272 C 15 Y(I) = 1.0
273 C Y(31) = 0.5
274 C T = 0.0
275 C TOUT = 0.1
276 C ISTATE = 0
277 C DO 30 IO = 1,4
278 C CALL DLSODIS (RESID, ADDASP, JACSP, NEQ, Y, YDOTI, T, TOUT,
279 C 1 ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RW, LRW, IW, LIW, MF)
280 C WRITE(6,20) T,IW(11),RW(11)
281 C 20 FORMAT(' At t =',F5.2,' No. steps =',I4,
282 C 1 ' Last step =',D12.4)
283 C IF (ISTATE .NE. 2) GO TO 90
284 C TOUT = TOUT + 0.1
285 C 30 CONTINUE
286 C WRITE (6,40) (Y(I),I=1,NEQ)
287 C 40 FORMAT(/' Final solution values..'/8(5D12.4/))
288 C WRITE(6,50) IW(17),IW(18),IW(11),IW(12),IW(13)
289 C NNZLU = IW(25) + IW(26) + NEQ
290 C WRITE(6,60) IW(19),NNZLU
291 C 50 FORMAT(/' Required RW size =',I5,' IW size =',I4/
292 C 1 ' No. steps =',I4,' No. r-s =',I4,' No. J-s =',i4)
293 C 60 FORMAT(' No. of nonzeros in P matrix =',I4,
294 C 1 ' No. of nonzeros in LU =',I4)
295 C STOP
296 C 90 WRITE (6,95) ISTATE
297 C 95 FORMAT(///' Error halt.. ISTATE =',I3)
298 C STOP
299 C END
300 C
301 C SUBROUTINE GFUN (N, T, Y, G)
302 C DOUBLE PRECISION T, Y, G, R4D, EODSQ
303 C DIMENSION G(N), Y(N)
304 C COMMON /TEST1/ R4D, EODSQ, NM1
305 C G(1) = R4D*(Y(N)**2-Y(2)**2) + EODSQ*(Y(2)-2.0*Y(1)+Y(N))
306 C DO 10 I = 2,NM1
307 C G(I) = R4D*(Y(I-1)**2 - Y(I+1)**2)
308 C 1 + EODSQ*(Y(I+1) - 2.0*Y(I) + Y(I-1))
309 C 10 CONTINUE
310 C G(N) = R4D*(Y(NM1)**2-Y(1)**2) + EODSQ*(Y(1)-2.0*Y(N)+Y(NM1))
311 C RETURN
312 C END
313 C
314 C SUBROUTINE RESID (N, T, Y, S, R, IRES)
315 C DOUBLE PRECISION T, Y, S, R, R4D, EODSQ
316 C DIMENSION Y(N), S(N), R(N)
317 C COMMON /TEST1/ R4D, EODSQ, NM1
318 C CALL GFUN (N, T, Y, R)
319 C R(1) = R(1) - (S(N) + 4.0*S(1) + S(2))/6.0
320 C DO 10 I = 2,NM1
321 C 10 R(I) = R(I) - (S(I-1) + 4.0*S(I) + S(I+1))/6.0
322 C R(N) = R(N) - (S(NM1) + 4.0*S(N) + S(1))/6.0
323 C RETURN
324 C END
325 C
326 C SUBROUTINE ADDASP (N, T, Y, J, IP, JP, P)
327 C DOUBLE PRECISION T, Y, P
328 C DIMENSION Y(N), IP(*), JP(*), P(N)
329 C JM1 = J - 1
330 C JP1 = J + 1
331 C IF (J .EQ. N) JP1 = 1
332 C IF (J .EQ. 1) JM1 = N
333 C P(J) = P(J) + (2.0/3.0)
334 C P(JP1) = P(JP1) + (1.0/6.0)
335 C P(JM1) = P(JM1) + (1.0/6.0)
336 C RETURN
337 C END
338 C
339 C SUBROUTINE JACSP (N, T, Y, S, J, IP, JP, PDJ)
340 C DOUBLE PRECISION T, Y, S, PDJ, R4D, EODSQ
341 C DIMENSION Y(N), S(N), IP(*), JP(*), PDJ(N)
342 C COMMON /TEST1/ R4D, EODSQ, NM1
343 C JM1 = J - 1
344 C JP1 = J + 1
345 C IF (J .EQ. 1) JM1 = N
346 C IF (J .EQ. N) JP1 = 1
347 C PDJ(JM1) = -2.0*R4D*Y(J) + EODSQ
348 C PDJ(J) = -2.0*EODSQ
349 C PDJ(JP1) = 2.0*R4D*Y(J) + EODSQ
350 C RETURN
351 C END
352 C
353 C The output of this program (on a CDC-7600 in single precision)
354 C is as follows:
355 C
356 C At t = 0.10 No. steps = 15 Last step = 1.6863e-02
357 C At t = 0.20 No. steps = 19 Last step = 2.4101e-02
358 C At t = 0.30 No. steps = 22 Last step = 4.3143e-02
359 C At t = 0.40 No. steps = 24 Last step = 5.7819e-02
360 C
361 C Final solution values..
362 C 1.8371e-02 1.3578e-02 1.5864e-02 2.3805e-02 3.7245e-02
363 C 5.6630e-02 8.2538e-02 1.1538e-01 1.5522e-01 2.0172e-01
364 C 2.5414e-01 3.1150e-01 3.7259e-01 4.3608e-01 5.0060e-01
365 C 5.6482e-01 6.2751e-01 6.8758e-01 7.4415e-01 7.9646e-01
366 C 8.4363e-01 8.8462e-01 9.1853e-01 9.4500e-01 9.6433e-01
367 C 9.7730e-01 9.8464e-01 9.8645e-01 9.8138e-01 9.6584e-01
368 C 9.3336e-01 8.7497e-01 7.8213e-01 6.5315e-01 4.9997e-01
369 C 3.4672e-01 2.1758e-01 1.2461e-01 6.6208e-02 3.3784e-02
370 C
371 C Required RW size = 1409 IW size = 30
372 C No. steps = 24 No. r-s = 33 No. J-s = 8
373 C No. of nonzeros in P matrix = 120 No. of nonzeros in LU = 194
374 C
375 C-----------------------------------------------------------------------
376 C Full Description of User Interface to DLSODIS.
377 C
378 C The user interface to DLSODIS consists of the following parts.
379 C
380 C 1. The call sequence to Subroutine DLSODIS, which is a driver
381 C routine for the solver. This includes descriptions of both
382 C the call sequence arguments and of user-supplied routines.
383 C Following these descriptions is a description of
384 C optional inputs available through the call sequence, and then
385 C a description of optional outputs (in the work arrays).
386 C
387 C 2. Descriptions of other routines in the DLSODIS package that may be
388 C (optionally) called by the user. These provide the ability to
389 C alter error message handling, save and restore the internal
390 C Common, and obtain specified derivatives of the solution y(t).
391 C
392 C 3. Descriptions of Common blocks to be declared in overlay
393 C or similar environments, or to be saved when doing an interrupt
394 C of the problem and continued solution later.
395 C
396 C 4. Description of two routines in the DLSODIS package, either of
397 C which the user may replace with his/her own version, if desired.
398 C These relate to the measurement of errors.
399 C
400 C-----------------------------------------------------------------------
401 C Part 1. Call Sequence.
402 C
403 C The call sequence parameters used for input only are
404 C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
405 C IOPT, LRW, LIW, MF,
406 C and those used for both input and output are
407 C Y, T, ISTATE, YDOTI.
408 C The work arrays RWORK and IWORK are also used for conditional and
409 C optional inputs and optional outputs. (The term output here refers
410 C to the return from Subroutine DLSODIS to the user's calling program.)
411 C
412 C The legality of input parameters will be thoroughly checked on the
413 C initial call for the problem, but not checked thereafter unless a
414 C change in input parameters is flagged by ISTATE = 3 on input.
415 C
416 C The descriptions of the call arguments are as follows.
417 C
418 C RES = the name of the user-supplied subroutine which supplies
419 C the residual vector for the ODE system, defined by
420 C r = g(t,y) - A(t,y) * s
421 C as a function of the scalar t and the vectors
422 C s and y (s approximates dy/dt). This subroutine
423 C is to have the form
424 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
425 C DOUBLE PRECISION T, Y(*), S(*), R(*)
426 C where NEQ, T, Y, S, and IRES are input, and R and
427 C IRES are output. Y, S, and R are arrays of length NEQ.
428 C On input, IRES indicates how DLSODIS will use the
429 C returned array R, as follows:
430 C IRES = 1 means that DLSODIS needs the full residual,
431 C r = g - A*s, exactly.
432 C IRES = -1 means that DLSODIS is using R only to compute
433 C the Jacobian dr/dy by difference quotients.
434 C The RES routine can ignore IRES, or it can omit some terms
435 C if IRES = -1. If A does not depend on y, then RES can
436 C just return R = g when IRES = -1. If g - A*s contains other
437 C additive terms that are independent of y, these can also be
438 C dropped, if done consistently, when IRES = -1.
439 C The subroutine should set the flag IRES if it
440 C encounters a halt condition or illegal input.
441 C Otherwise, it should not reset IRES. On output,
442 C IRES = 1 or -1 represents a normal return, and
443 C DLSODIS continues integrating the ODE. Leave IRES
444 C unchanged from its input value.
445 C IRES = 2 tells DLSODIS to immediately return control
446 C to the calling program, with ISTATE = 3. This lets
447 C the calling program change parameters of the problem
448 C if necessary.
449 C IRES = 3 represents an error condition (for example, an
450 C illegal value of y). DLSODIS tries to integrate the system
451 C without getting IRES = 3 from RES. If it cannot, DLSODIS
452 C returns with ISTATE = -7 or -1.
453 C On a return with ISTATE = 3, -1, or -7, the values
454 C of T and Y returned correspond to the last point reached
455 C successfully without getting the flag IRES = 2 or 3.
456 C The flag values IRES = 2 and 3 should not be used to
457 C handle switches or root-stop conditions. This is better
458 C done by calling DLSODIS in a one-step mode and checking the
459 C stopping function for a sign change at each step.
460 C If quantities computed in the RES routine are needed
461 C externally to DLSODIS, an extra call to RES should be made
462 C for this purpose, for consistent and accurate results.
463 C To get the current dy/dt for the S argument, use DINTDY.
464 C RES must be declared External in the calling
465 C program. See note below for more about RES.
466 C
467 C ADDA = the name of the user-supplied subroutine which adds the
468 C matrix A = A(t,y) to another matrix stored in sparse form.
469 C This subroutine is to have the form
470 C SUBROUTINE ADDA (NEQ, T, Y, J, IAN, JAN, P)
471 C DOUBLE PRECISION T, Y(*), P(*)
472 C INTEGER IAN(*), JAN(*)
473 C where NEQ, T, Y, J, IAN, JAN, and P are input. This routine
474 C should add the J-th column of matrix A to the array P, of
475 C length NEQ. Thus a(i,J) is to be added to P(i) for all
476 C relevant values of i. Here T and Y have the same meaning as
477 C in Subroutine RES, and J is a column index (1 to NEQ).
478 C IAN and JAN are undefined in calls to ADDA for structure
479 C determination (MOSS .ne. 0). Otherwise, IAN and JAN are
480 C structure descriptors, as defined under optional outputs
481 C below, and so can be used to determine the relevant row
482 C indices i, if desired.
483 C Calls to ADDA are made with J = 1,...,NEQ, in that
484 C order. ADDA must not alter its input arguments.
485 C ADDA must be declared External in the calling program.
486 C See note below for more information about ADDA.
487 C
488 C JAC = the name of the user-supplied subroutine which supplies
489 C the Jacobian matrix, dr/dy, where r = g - A*s. JAC is
490 C required if MITER = 1, or MOSS = 1 or 3. Otherwise a dummy
491 C name can be passed. This subroutine is to have the form
492 C SUBROUTINE JAC (NEQ, T, Y, S, J, IAN, JAN, PDJ)
493 C DOUBLE PRECISION T, Y(*), S(*), PDJ(*)
494 C INTEGER IAN(*), JAN(*)
495 C where NEQ, T, Y, S, J, IAN, and JAN are input. The
496 C array PDJ, of length NEQ, is to be loaded with column J
497 C of the Jacobian on output. Thus dr(i)/dy(J) is to be
498 C loaded into PDJ(i) for all relevant values of i.
499 C Here T, Y, and S have the same meaning as in Subroutine RES,
500 C and J is a column index (1 to NEQ). IAN and JAN
501 C are undefined in calls to JAC for structure determination
502 C (MOSS .ne. 0). Otherwise, IAN and JAN are structure
503 C descriptors, as defined under optional outputs below, and
504 C so can be used to determine the relevant row indices i, if
505 C desired.
506 C JAC need not provide dr/dy exactly. A crude
507 C approximation (possibly with greater sparsity) will do.
508 C In any case, PDJ is preset to zero by the solver,
509 C so that only the nonzero elements need be loaded by JAC.
510 C Calls to JAC are made with J = 1,...,NEQ, in that order, and
511 C each such set of calls is preceded by a call to RES with the
512 C same arguments NEQ, T, Y, S, and IRES. Thus to gain some
513 C efficiency intermediate quantities shared by both calculations
514 C may be saved in a user Common block by RES and not recomputed
515 C by JAC, if desired. JAC must not alter its input arguments.
516 C JAC must be declared External in the calling program.
517 C See note below for more about JAC.
518 C
519 C Note on RES, ADDA, and JAC:
520 C These subroutines may access user-defined quantities in
521 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
522 C (dimensioned in the subroutines) and/or Y has length
523 C exceeding NEQ(1). However, these subroutines should not
524 C alter NEQ(1), Y(1),...,Y(NEQ) or any other input variables.
525 C See the descriptions of NEQ and Y below.
526 C
527 C NEQ = the size of the system (number of first order ordinary
528 C differential equations or scalar algebraic equations).
529 C Used only for input.
530 C NEQ may be decreased, but not increased, during the problem.
531 C If NEQ is decreased (with ISTATE = 3 on input), the
532 C remaining components of Y should be left undisturbed, if
533 C these are to be accessed in RES, ADDA, or JAC.
534 C
535 C Normally, NEQ is a scalar, and it is generally referred to
536 C as a scalar in this user interface description. However,
537 C NEQ may be an array, with NEQ(1) set to the system size.
538 C (The DLSODIS package accesses only NEQ(1).) In either case,
539 C this parameter is passed as the NEQ argument in all calls
540 C to RES, ADDA, and JAC. Hence, if it is an array,
541 C locations NEQ(2),... may be used to store other integer data
542 C and pass it to RES, ADDA, or JAC. Each such subroutine
543 C must include NEQ in a Dimension statement in that case.
544 C
545 C Y = a real array for the vector of dependent variables, of
546 C length NEQ or more. Used for both input and output on the
547 C first call (ISTATE = 0 or 1), and only for output on other
548 C calls. On the first call, Y must contain the vector of
549 C initial values. On output, Y contains the computed solution
550 C vector, evaluated at T. If desired, the Y array may be used
551 C for other purposes between calls to the solver.
552 C
553 C This array is passed as the Y argument in all calls to RES,
554 C ADDA, and JAC. Hence its length may exceed NEQ,
555 C and locations Y(NEQ+1),... may be used to store other real
556 C data and pass it to RES, ADDA, or JAC. (The DLSODIS
557 C package accesses only Y(1),...,Y(NEQ). )
558 C
559 C YDOTI = a real array for the initial value of the vector
560 C dy/dt and for work space, of dimension at least NEQ.
561 C
562 C On input:
563 C If ISTATE = 0 then DLSODIS will compute the initial value
564 C of dy/dt, if A is nonsingular. Thus YDOTI will
565 C serve only as work space and may have any value.
566 C If ISTATE = 1 then YDOTI must contain the initial value
567 C of dy/dt.
568 C If ISTATE = 2 or 3 (continuation calls) then YDOTI
569 C may have any value.
570 C Note: If the initial value of A is singular, then
571 C DLSODIS cannot compute the initial value of dy/dt, so
572 C it must be provided in YDOTI, with ISTATE = 1.
573 C
574 C On output, when DLSODIS terminates abnormally with ISTATE =
575 C -1, -4, or -5, YDOTI will contain the residual
576 C r = g(t,y) - A(t,y)*(dy/dt). If r is large, t is near
577 C its initial value, and YDOTI is supplied with ISTATE = 1,
578 C there may have been an incorrect input value of
579 C YDOTI = dy/dt, or the problem (as given to DLSODIS)
580 C may not have a solution.
581 C
582 C If desired, the YDOTI array may be used for other
583 C purposes between calls to the solver.
584 C
585 C T = the independent variable. On input, T is used only on the
586 C first call, as the initial point of the integration.
587 C On output, after each call, T is the value at which a
588 C computed solution y is evaluated (usually the same as TOUT).
589 C On an error return, T is the farthest point reached.
590 C
591 C TOUT = the next value of t at which a computed solution is desired.
592 C Used only for input.
593 C
594 C When starting the problem (ISTATE = 0 or 1), TOUT may be
595 C equal to T for one call, then should .ne. T for the next
596 C call. For the initial T, an input value of TOUT .ne. T is
597 C used in order to determine the direction of the integration
598 C (i.e. the algebraic sign of the step sizes) and the rough
599 C scale of the problem. Integration in either direction
600 C (forward or backward in t) is permitted.
601 C
602 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
603 C the first call (i.e. the first call with TOUT .ne. T).
604 C Otherwise, TOUT is required on every call.
605 C
606 C If ITASK = 1, 3, or 4, the values of TOUT need not be
607 C monotone, but a value of TOUT which backs up is limited
608 C to the current internal T interval, whose endpoints are
609 C TCUR - HU and TCUR (see optional outputs, below, for
610 C TCUR and HU).
611 C
612 C ITOL = an indicator for the type of error control. See
613 C description below under ATOL. Used only for input.
614 C
615 C RTOL = a relative error tolerance parameter, either a scalar or
616 C an array of length NEQ. See description below under ATOL.
617 C Input only.
618 C
619 C ATOL = an absolute error tolerance parameter, either a scalar or
620 C an array of length NEQ. Input only.
621 C
622 C The input parameters ITOL, RTOL, and ATOL determine
623 C the error control performed by the solver. The solver will
624 C control the vector E = (E(i)) of estimated local errors
625 C in y, according to an inequality of the form
626 C RMS-norm of ( E(i)/EWT(i) ) .le. 1,
627 C where EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
628 C and the RMS-norm (root-mean-square norm) here is
629 C RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
630 C is a vector of weights which must always be positive, and
631 C the values of RTOL and ATOL should all be non-negative.
632 C The following table gives the types (scalar/array) of
633 C RTOL and ATOL, and the corresponding form of EWT(i).
634 C
635 C ITOL RTOL ATOL EWT(i)
636 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
637 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
638 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
639 C 4 array scalar RTOL(i)*ABS(Y(i)) + ATOL(i)
640 C
641 C When either of these parameters is a scalar, it need not
642 C be dimensioned in the user's calling program.
643 C
644 C If none of the above choices (with ITOL, RTOL, and ATOL
645 C fixed throughout the problem) is suitable, more general
646 C error controls can be obtained by substituting
647 C user-supplied routines for the setting of EWT and/or for
648 C the norm calculation. See Part 4 below.
649 C
650 C If global errors are to be estimated by making a repeated
651 C run on the same problem with smaller tolerances, then all
652 C components of RTOL and ATOL (i.e. of EWT) should be scaled
653 C down uniformly.
654 C
655 C ITASK = an index specifying the task to be performed.
656 C Input only. ITASK has the following values and meanings.
657 C 1 means normal computation of output values of y(t) at
658 C t = TOUT (by overshooting and interpolating).
659 C 2 means take one step only and return.
660 C 3 means stop at the first internal mesh point at or
661 C beyond t = TOUT and return.
662 C 4 means normal computation of output values of y(t) at
663 C t = TOUT but without overshooting t = TCRIT.
664 C TCRIT must be input as RWORK(1). TCRIT may be equal to
665 C or beyond TOUT, but not behind it in the direction of
666 C integration. This option is useful if the problem
667 C has a singularity at or beyond t = TCRIT.
668 C 5 means take one step, without passing TCRIT, and return.
669 C TCRIT must be input as RWORK(1).
670 C
671 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
672 C (within roundoff), it will return T = TCRIT (exactly) to
673 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
674 C in which case answers at t = TOUT are returned first).
675 C
676 C ISTATE = an index used for input and output to specify the
677 C state of the calculation.
678 C
679 C On input, the values of ISTATE are as follows.
680 C 0 means this is the first call for the problem, and
681 C DLSODIS is to compute the initial value of dy/dt
682 C (while doing other initializations). See note below.
683 C 1 means this is the first call for the problem, and
684 C the initial value of dy/dt has been supplied in
685 C YDOTI (DLSODIS will do other initializations).
686 C See note below.
687 C 2 means this is not the first call, and the calculation
688 C is to continue normally, with no change in any input
689 C parameters except possibly TOUT and ITASK.
690 C (If ITOL, RTOL, and/or ATOL are changed between calls
691 C with ISTATE = 2, the new values will be used but not
692 C tested for legality.)
693 C 3 means this is not the first call, and the
694 C calculation is to continue normally, but with
695 C a change in input parameters other than
696 C TOUT and ITASK. Changes are allowed in
697 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF,
698 C the conditional inputs IA, JA, IC, and JC,
699 C and any of the optional inputs except H0.
700 C A call with ISTATE = 3 will cause the sparsity
701 C structure of the problem to be recomputed.
702 C (Structure information is reread from IA and JA if
703 C MOSS = 0, 3, or 4 and from IC and JC if MOSS = 0).
704 C Note: A preliminary call with TOUT = T is not counted
705 C as a first call here, as no initialization or checking of
706 C input is done. (Such a call is sometimes useful for the
707 C purpose of outputting the initial conditions.)
708 C Thus the first call for which TOUT .ne. T requires
709 C ISTATE = 0 or 1 on input.
710 C
711 C On output, ISTATE has the following values and meanings.
712 C 0 or 1 means nothing was done; TOUT = T and
713 C ISTATE = 0 or 1 on input.
714 C 2 means that the integration was performed successfully.
715 C 3 means that the user-supplied Subroutine RES signalled
716 C DLSODIS to halt the integration and return (IRES = 2).
717 C Integration as far as T was achieved with no occurrence
718 C of IRES = 2, but this flag was set on attempting the
719 C next step.
720 C -1 means an excessive amount of work (more than MXSTEP
721 C steps) was done on this call, before completing the
722 C requested task, but the integration was otherwise
723 C successful as far as T. (MXSTEP is an optional input
724 C and is normally 500.) To continue, the user may
725 C simply reset ISTATE to a value .gt. 1 and call again
726 C (the excess work step counter will be reset to 0).
727 C In addition, the user may increase MXSTEP to avoid
728 C this error return (see below on optional inputs).
729 C -2 means too much accuracy was requested for the precision
730 C of the machine being used. This was detected before
731 C completing the requested task, but the integration
732 C was successful as far as T. To continue, the tolerance
733 C parameters must be reset, and ISTATE must be set
734 C to 3. The optional output TOLSF may be used for this
735 C purpose. (Note: If this condition is detected before
736 C taking any steps, then an illegal input return
737 C (ISTATE = -3) occurs instead.)
738 C -3 means illegal input was detected, before taking any
739 C integration steps. See written message for details.
740 C Note: If the solver detects an infinite loop of calls
741 C to the solver with illegal input, it will cause
742 C the run to stop.
743 C -4 means there were repeated error test failures on
744 C one attempted step, before completing the requested
745 C task, but the integration was successful as far as T.
746 C The problem may have a singularity, or the input
747 C may be inappropriate.
748 C -5 means there were repeated convergence test failures on
749 C one attempted step, before completing the requested
750 C task, but the integration was successful as far as T.
751 C This may be caused by an inaccurate Jacobian matrix.
752 C -6 means EWT(i) became zero for some i during the
753 C integration. Pure relative error control (ATOL(i) = 0.0)
754 C was requested on a variable which has now vanished.
755 C the integration was successful as far as T.
756 C -7 means that the user-supplied Subroutine RES set
757 C its error flag (IRES = 3) despite repeated tries by
758 C DLSODIS to avoid that condition.
759 C -8 means that ISTATE was 0 on input but DLSODIS was unable
760 C to compute the initial value of dy/dt. See the
761 C printed message for details.
762 C -9 means a fatal error return flag came from the sparse
763 C solver CDRV by way of DPRJIS or DSOLSS (numerical
764 C factorization or backsolve). This should never happen.
765 C The integration was successful as far as T.
766 C
767 C Note: An error return with ISTATE = -1, -4, or -5
768 C may mean that the sparsity structure of the
769 C problem has changed significantly since it was last
770 C determined (or input). In that case, one can attempt to
771 C complete the integration by setting ISTATE = 3 on the next
772 C call, so that a new structure determination is done.
773 C
774 C Note: Since the normal output value of ISTATE is 2,
775 C it does not need to be reset for normal continuation.
776 C similarly, ISTATE (= 3) need not be reset if RES told
777 C DLSODIS to return because the calling program must change
778 C the parameters of the problem.
779 C Also, since a negative input value of ISTATE will be
780 C regarded as illegal, a negative output value requires the
781 C user to change it, and possibly other inputs, before
782 C calling the solver again.
783 C
784 C IOPT = an integer flag to specify whether or not any optional
785 C inputs are being used on this call. Input only.
786 C The optional inputs are listed separately below.
787 C IOPT = 0 means no optional inputs are being used.
788 C Default values will be used in all cases.
789 C IOPT = 1 means one or more optional inputs are being used.
790 C
791 C RWORK = a work array used for a mixture of real (double precision)
792 C and integer work space.
793 C The length of RWORK (in real words) must be at least
794 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM where
795 C NYH = the initial value of NEQ,
796 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
797 C smaller value is given as an optional input),
798 C LWM = 2*NNZ + 2*NEQ + (NNZ+9*NEQ)/LENRAT if MITER = 1,
799 C LWM = 2*NNZ + 2*NEQ + (NNZ+10*NEQ)/LENRAT if MITER = 2.
800 C in the above formulas,
801 C NNZ = number of nonzero elements in the iteration matrix
802 C P = A - con*J (con is a constant and J is the
803 C Jacobian matrix dr/dy).
804 C LENRAT = the real to integer wordlength ratio (usually 1 in
805 C single precision and 2 in double precision).
806 C (See the MF description for METH and MITER.)
807 C Thus if MAXORD has its default value and NEQ is constant,
808 C the minimum length of RWORK is:
809 C 20 + 16*NEQ + LWM for MF = 11, 111, 311, 12, 212, 412,
810 C 20 + 9*NEQ + LWM for MF = 21, 121, 321, 22, 222, 422.
811 C The above formula for LWM is only a crude lower bound.
812 C The required length of RWORK cannot be readily predicted
813 C in general, as it depends on the sparsity structure
814 C of the problem. Some experimentation may be necessary.
815 C
816 C The first 20 words of RWORK are reserved for conditional
817 C and optional inputs and optional outputs.
818 C
819 C The following word in RWORK is a conditional input:
820 C RWORK(1) = TCRIT = critical value of t which the solver
821 C is not to overshoot. Required if ITASK is
822 C 4 or 5, and ignored otherwise. (See ITASK.)
823 C
824 C LRW = the length of the array RWORK, as declared by the user.
825 C (This will be checked by the solver.)
826 C
827 C IWORK = an integer work array. The length of IWORK must be at least
828 C 32 + 2*NEQ + NZA + NZC for MOSS = 0,
829 C 30 for MOSS = 1 or 2,
830 C 31 + NEQ + NZA for MOSS = 3 or 4.
831 C (NZA is the number of nonzero elements in matrix A, and
832 C NZC is the number of nonzero elements in dr/dy.)
833 C
834 C In DLSODIS, IWORK is used for conditional and
835 C optional inputs and optional outputs.
836 C
837 C The following two blocks of words in IWORK are conditional
838 C inputs, required if MOSS = 0, 3, or 4, but not otherwise
839 C (see the description of MF for MOSS).
840 C IWORK(30+j) = IA(j) (j=1,...,NEQ+1)
841 C IWORK(31+NEQ+k) = JA(k) (k=1,...,NZA)
842 C The two arrays IA and JA describe the sparsity structure
843 C to be assumed for the matrix A. JA contains the row
844 C indices where nonzero elements occur, reading in columnwise
845 C order, and IA contains the starting locations in JA of the
846 C descriptions of columns 1,...,NEQ, in that order, with
847 C IA(1) = 1. Thus, for each column index j = 1,...,NEQ, the
848 C values of the row index i in column j where a nonzero
849 C element may occur are given by
850 C i = JA(k), where IA(j) .le. k .lt. IA(j+1).
851 C If NZA is the total number of nonzero locations assumed,
852 C then the length of the JA array is NZA, and IA(NEQ+1) must
853 C be NZA + 1. Duplicate entries are not allowed.
854 C The following additional blocks of words are required
855 C if MOSS = 0, but not otherwise. If LC = 31 + NEQ + NZA, then
856 C IWORK(LC+j) = IC(j) (j=1,...,NEQ+1), and
857 C IWORK(LC+NEQ+1+k) = JC(k) (k=1,...,NZC)
858 C The two arrays IC and JC describe the sparsity
859 C structure to be assumed for the Jacobian matrix dr/dy.
860 C They are used in the same manner as the above IA and JA
861 C arrays. If NZC is the number of nonzero locations
862 C assumed, then the length of the JC array is NZC, and
863 C IC(NEQ+1) must be NZC + 1. Duplicate entries are not
864 C allowed.
865 C
866 C LIW = the length of the array IWORK, as declared by the user.
867 C (This will be checked by the solver.)
868 C
869 C Note: The work arrays must not be altered between calls to DLSODIS
870 C for the same problem, except possibly for the conditional and
871 C optional inputs, and except for the last 3*NEQ words of RWORK.
872 C The latter space is used for internal scratch space, and so is
873 C available for use by the user outside DLSODIS between calls, if
874 C desired (but not for use by RES, ADDA, or JAC).
875 C
876 C MF = the method flag. Used only for input.
877 C MF has three decimal digits-- MOSS, METH, and MITER.
878 C For standard options:
879 C MF = 100*MOSS + 10*METH + MITER.
880 C MOSS indicates the method to be used to obtain the sparsity
881 C structure of the Jacobian matrix:
882 C MOSS = 0 means the user has supplied IA, JA, IC, and JC
883 C (see descriptions under IWORK above).
884 C MOSS = 1 means the user has supplied JAC (see below) and
885 C the structure will be obtained from NEQ initial
886 C calls to JAC and NEQ initial calls to ADDA.
887 C MOSS = 2 means the structure will be obtained from NEQ+1
888 C initial calls to RES and NEQ initial calls to ADDA
889 C MOSS = 3 like MOSS = 1, except user has supplied IA and JA.
890 C MOSS = 4 like MOSS = 2, except user has supplied IA and JA.
891 C METH indicates the basic linear multistep method:
892 C METH = 1 means the implicit Adams method.
893 C METH = 2 means the method based on Backward
894 C Differentiation Formulas (BDFs).
895 C The BDF method is strongly preferred for stiff problems,
896 C while the Adams method is preferred when the problem is
897 C not stiff. If the matrix A(t,y) is nonsingular,
898 C stiffness here can be taken to mean that of the explicit
899 C ODE system dy/dt = A-inverse * g. If A is singular,
900 C the concept of stiffness is not well defined.
901 C If you do not know whether the problem is stiff, we
902 C recommend using METH = 2. If it is stiff, the advantage
903 C of METH = 2 over METH = 1 will be great, while if it is
904 C not stiff, the advantage of METH = 1 will be slight.
905 C If maximum efficiency is important, some experimentation
906 C with METH may be necessary.
907 C MITER indicates the corrector iteration method:
908 C MITER = 1 means chord iteration with a user-supplied
909 C sparse Jacobian, given by Subroutine JAC.
910 C MITER = 2 means chord iteration with an internally
911 C generated (difference quotient) sparse
912 C Jacobian (using NGP extra calls to RES per
913 C dr/dy value, where NGP is an optional
914 C output described below.)
915 C If MITER = 1 or MOSS = 1 or 3 the user must supply a
916 C Subroutine JAC (the name is arbitrary) as described above
917 C under JAC. Otherwise, a dummy argument can be used.
918 C
919 C The standard choices for MF are:
920 C MF = 21 or 22 for a stiff problem with IA/JA and IC/JC
921 C supplied,
922 C MF = 121 for a stiff problem with JAC supplied, but not
923 C IA/JA or IC/JC,
924 C MF = 222 for a stiff problem with neither IA/JA, IC/JC/,
925 C nor JAC supplied,
926 C MF = 321 for a stiff problem with IA/JA and JAC supplied,
927 C but not IC/JC,
928 C MF = 422 for a stiff problem with IA/JA supplied, but not
929 C IC/JC or JAC.
930 C
931 C The sparseness structure can be changed during the problem
932 C by making a call to DLSODIS with ISTATE = 3.
933 C-----------------------------------------------------------------------
934 C Optional Inputs.
935 C
936 C The following is a list of the optional inputs provided for in the
937 C call sequence. (See also Part 2.) For each such input variable,
938 C this table lists its name as used in this documentation, its
939 C location in the call sequence, its meaning, and the default value.
940 C The use of any of these inputs requires IOPT = 1, and in that
941 C case all of these inputs are examined. A value of zero for any
942 C of these optional inputs will cause the default value to be used.
943 C Thus to use a subset of the optional inputs, simply preload
944 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
945 C then set those of interest to nonzero values.
946 C
947 C Name Location Meaning and Default Value
948 C
949 C H0 RWORK(5) the step size to be attempted on the first step.
950 C The default value is determined by the solver.
951 C
952 C HMAX RWORK(6) the maximum absolute step size allowed.
953 C The default value is infinite.
954 C
955 C HMIN RWORK(7) the minimum absolute step size allowed.
956 C The default value is 0. (This lower bound is not
957 C enforced on the final step before reaching TCRIT
958 C when ITASK = 4 or 5.)
959 C
960 C MAXORD IWORK(5) the maximum order to be allowed. The default
961 C value is 12 if METH = 1, and 5 if METH = 2.
962 C If MAXORD exceeds the default value, it will
963 C be reduced to the default value.
964 C If MAXORD is changed during the problem, it may
965 C cause the current order to be reduced.
966 C
967 C MXSTEP IWORK(6) maximum number of (internally defined) steps
968 C allowed during one call to the solver.
969 C The default value is 500.
970 C
971 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
972 C warning that T + H = T on a step (H = step size).
973 C This must be positive to result in a non-default
974 C value. The default value is 10.
975 C-----------------------------------------------------------------------
976 C Optional Outputs.
977 C
978 C As optional additional output from DLSODIS, the variables listed
979 C below are quantities related to the performance of DLSODIS
980 C which are available to the user. These are communicated by way of
981 C the work arrays, but also have internal mnemonic names as shown.
982 C Except where stated otherwise, all of these outputs are defined
983 C on any successful return from DLSODIS, and on any return with
984 C ISTATE = -1, -2, -4, -5, -6, or -7. On a return with -3 (illegal
985 C input) or -8, they will be unchanged from their existing values
986 C (if any), except possibly for TOLSF, LENRW, and LENIW.
987 C On any error return, outputs relevant to the error will be defined,
988 C as noted below.
989 C
990 C Name Location Meaning
991 C
992 C HU RWORK(11) the step size in t last used (successfully).
993 C
994 C HCUR RWORK(12) the step size to be attempted on the next step.
995 C
996 C TCUR RWORK(13) the current value of the independent variable
997 C which the solver has actually reached, i.e. the
998 C current internal mesh point in t. On output, TCUR
999 C will always be at least as far as the argument
1000 C T, but may be farther (if interpolation was done).
1001 C
1002 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
1003 C computed when a request for too much accuracy was
1004 C detected (ISTATE = -3 if detected at the start of
1005 C the problem, ISTATE = -2 otherwise). If ITOL is
1006 C left unaltered but RTOL and ATOL are uniformly
1007 C scaled up by a factor of TOLSF for the next call,
1008 C then the solver is deemed likely to succeed.
1009 C (The user may also ignore TOLSF and alter the
1010 C tolerance parameters in any other way appropriate.)
1011 C
1012 C NST IWORK(11) the number of steps taken for the problem so far.
1013 C
1014 C NRE IWORK(12) the number of residual evaluations (RES calls)
1015 C for the problem so far, excluding those for
1016 C structure determination (MOSS = 2 or 4).
1017 C
1018 C NJE IWORK(13) the number of Jacobian evaluations (each involving
1019 C an evaluation of A and dr/dy) for the problem so
1020 C far, excluding those for structure determination
1021 C (MOSS = 1 or 3). This equals the number of calls
1022 C to ADDA and (if MITER = 1) JAC.
1023 C
1024 C NQU IWORK(14) the method order last used (successfully).
1025 C
1026 C NQCUR IWORK(15) the order to be attempted on the next step.
1027 C
1028 C IMXER IWORK(16) the index of the component of largest magnitude in
1029 C the weighted local error vector ( E(i)/EWT(i) ),
1030 C on an error return with ISTATE = -4 or -5.
1031 C
1032 C LENRW IWORK(17) the length of RWORK actually required.
1033 C This is defined on normal returns and on an illegal
1034 C input return for insufficient storage.
1035 C
1036 C LENIW IWORK(18) the length of IWORK actually required.
1037 C This is defined on normal returns and on an illegal
1038 C input return for insufficient storage.
1039 C
1040 C NNZ IWORK(19) the number of nonzero elements in the iteration
1041 C matrix P = A - con*J (con is a constant and
1042 C J is the Jacobian matrix dr/dy).
1043 C
1044 C NGP IWORK(20) the number of groups of column indices, used in
1045 C difference quotient Jacobian aproximations if
1046 C MITER = 2. This is also the number of extra RES
1047 C evaluations needed for each Jacobian evaluation.
1048 C
1049 C NLU IWORK(21) the number of sparse LU decompositions for the
1050 C problem so far. (Excludes the LU decomposition
1051 C necessary when ISTATE = 0.)
1052 C
1053 C LYH IWORK(22) the base address in RWORK of the history array YH,
1054 C described below in this list.
1055 C
1056 C IPIAN IWORK(23) the base address of the structure descriptor array
1057 C IAN, described below in this list.
1058 C
1059 C IPJAN IWORK(24) the base address of the structure descriptor array
1060 C JAN, described below in this list.
1061 C
1062 C NZL IWORK(25) the number of nonzero elements in the strict lower
1063 C triangle of the LU factorization used in the chord
1064 C iteration.
1065 C
1066 C NZU IWORK(26) the number of nonzero elements in the strict upper
1067 C triangle of the LU factorization used in the chord
1068 C iteration. The total number of nonzeros in the
1069 C factorization is therefore NZL + NZU + NEQ.
1070 C
1071 C The following four arrays are segments of the RWORK array which
1072 C may also be of interest to the user as optional outputs.
1073 C For each array, the table below gives its internal name,
1074 C its base address, and its description.
1075 C For YH and ACOR, the base addresses are in RWORK (a real array).
1076 C The integer arrays IAN and JAN are to be obtained by declaring an
1077 C integer array IWK and identifying IWK(1) with RWORK(21), using either
1078 C an equivalence statement or a subroutine call. Then the base
1079 C addresses IPIAN (of IAN) and IPJAN (of JAN) in IWK are to be obtained
1080 C as optional outputs IWORK(23) and IWORK(24), respectively.
1081 C Thus IAN(1) is IWK(ipian), etc.
1082 C
1083 C Name Base Address Description
1084 C
1085 C IAN IPIAN (in IWK) structure descriptor array of size NEQ + 1.
1086 C JAN IPJAN (in IWK) structure descriptor array of size NNZ.
1087 C (see above) IAN and JAN together describe the sparsity
1088 C structure of the iteration matrix
1089 C P = A - con*J, as used by DLSODIS.
1090 C JAN contains the row indices of the nonzero
1091 C locations, reading in columnwise order, and
1092 C IAN contains the starting locations in JAN of
1093 C the descriptions of columns 1,...,NEQ, in
1094 C that order, with IAN(1) = 1. Thus for each
1095 C j = 1,...,NEQ, the row indices i of the
1096 C nonzero locations in column j are
1097 C i = JAN(k), IAN(j) .le. k .lt. IAN(j+1).
1098 C Note that IAN(NEQ+1) = NNZ + 1.
1099 C YH LYH the Nordsieck history array, of size NYH by
1100 C (optional (NQCUR + 1), where NYH is the initial value
1101 C output) of NEQ. For j = 0,1,...,NQCUR, column j+1
1102 C of YH contains HCUR**j/factorial(j) times
1103 C the j-th derivative of the interpolating
1104 C polynomial currently representing the solution,
1105 C evaluated at t = TCUR. The base address LYH
1106 C is another optional output, listed above.
1107 C
1108 C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
1109 C corrections on each step, scaled on output to
1110 C represent the estimated local error in y on the
1111 C last step. This is the vector E in the
1112 C description of the error control. It is defined
1113 C only on a return from DLSODIS with ISTATE = 2.
1114 C
1115 C-----------------------------------------------------------------------
1116 C Part 2. Other Routines Callable.
1117 C
1118 C The following are optional calls which the user may make to
1119 C gain additional capabilities in conjunction with DLSODIS.
1120 C (The routines XSETUN and XSETF are designed to conform to the
1121 C SLATEC error handling package.)
1122 C
1123 C Form of Call Function
1124 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
1125 C output of messages from DLSODIS, if
1126 C The default is not desired.
1127 C The default value of LUN is 6.
1128 C
1129 C CALL XSETF(MFLAG) Set a flag to control the printing of
1130 C messages by DLSODIS.
1131 C MFLAG = 0 means do not print. (Danger:
1132 C This risks losing valuable information.)
1133 C MFLAG = 1 means print (the default).
1134 C
1135 C Either of the above calls may be made at
1136 C any time and will take effect immediately.
1137 C
1138 C CALL DSRCMS(RSAV,ISAV,JOB) saves and restores the contents of
1139 C the internal Common blocks used by
1140 C DLSODIS (see Part 3 below).
1141 C RSAV must be a real array of length 224
1142 C or more, and ISAV must be an integer
1143 C array of length 71 or more.
1144 C JOB=1 means save Common into RSAV/ISAV.
1145 C JOB=2 means restore Common from RSAV/ISAV.
1146 C DSRCMS is useful if one is
1147 C interrupting a run and restarting
1148 C later, or alternating between two or
1149 C more problems solved with DLSODIS.
1150 C
1151 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
1152 C (see below) orders, at a specified point t, if
1153 C desired. It may be called only after
1154 C a successful return from DLSODIS.
1155 C
1156 C The detailed instructions for using DINTDY are as follows.
1157 C The form of the call is:
1158 C
1159 C LYH = IWORK(22)
1160 C CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG)
1161 C
1162 C The input parameters are:
1163 C
1164 C T = value of independent variable where answers are desired
1165 C (normally the same as the T last returned by DLSODIS).
1166 C For valid results, T must lie between TCUR - HU and TCUR.
1167 C (See optional outputs for TCUR and HU.)
1168 C K = integer order of the derivative desired. K must satisfy
1169 C 0 .le. K .le. NQCUR, where NQCUR is the current order
1170 C (see optional outputs). The capability corresponding
1171 C to K = 0, i.e. computing y(t), is already provided
1172 C by DLSODIS directly. Since NQCUR .ge. 1, the first
1173 C derivative dy/dt is always available with DINTDY.
1174 C LYH = the base address of the history array YH, obtained
1175 C as an optional output as shown above.
1176 C NYH = column length of YH, equal to the initial value of NEQ.
1177 C
1178 C The output parameters are:
1179 C
1180 C DKY = a real array of length NEQ containing the computed value
1181 C of the K-th derivative of y(t).
1182 C IFLAG = integer flag, returned as 0 if K and T were legal,
1183 C -1 if K was illegal, and -2 if T was illegal.
1184 C On an error return, a message is also written.
1185 C-----------------------------------------------------------------------
1186 C Part 3. Common Blocks.
1187 C
1188 C If DLSODIS is to be used in an overlay situation, the user
1189 C must declare, in the primary overlay, the variables in:
1190 C (1) the call sequence to DLSODIS, and
1191 C (2) the two internal Common blocks
1192 C /DLS001/ of length 255 (218 double precision words
1193 C followed by 37 integer words),
1194 C /DLSS01/ of length 40 (6 double precision words
1195 C followed by 34 integer words).
1196 C
1197 C If DLSODIS is used on a system in which the contents of internal
1198 C Common blocks are not preserved between calls, the user should
1199 C declare the above Common blocks in the calling program to insure
1200 C that their contents are preserved.
1201 C
1202 C If the solution of a given problem by DLSODIS is to be interrupted
1203 C and then later continued, such as when restarting an interrupted run
1204 C or alternating between two or more problems, the user should save,
1205 C following the return from the last DLSODIS call prior to the
1206 C interruption, the contents of the call sequence variables and the
1207 C internal Common blocks, and later restore these values before the
1208 C next DLSODIS call for that problem. To save and restore the Common
1209 C blocks, use Subroutines DSRCMS (see Part 2 above).
1210 C
1211 C-----------------------------------------------------------------------
1212 C Part 4. Optionally Replaceable Solver Routines.
1213 C
1214 C Below are descriptions of two routines in the DLSODIS package which
1215 C relate to the measurement of errors. Either routine can be
1216 C replaced by a user-supplied version, if desired. However, since such
1217 C a replacement may have a major impact on performance, it should be
1218 C done only when absolutely necessary, and only with great caution.
1219 C (Note: The means by which the package version of a routine is
1220 C superseded by the user's version may be system-dependent.)
1221 C
1222 C (a) DEWSET.
1223 C The following subroutine is called just before each internal
1224 C integration step, and sets the array of error weights, EWT, as
1225 C described under ITOL/RTOL/ATOL above:
1226 C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1227 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSODIS call sequence,
1228 C YCUR contains the current dependent variable vector, and
1229 C EWT is the array of weights set by DEWSET.
1230 C
1231 C If the user supplies this subroutine, it must return in EWT(i)
1232 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1233 C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
1234 C routine (see below), and also used by DLSODIS in the computation
1235 C of the optional output IMXER, and the increments for difference
1236 C quotient Jacobians.
1237 C
1238 C In the user-supplied version of DEWSET, it may be desirable to use
1239 C the current values of derivatives of y. Derivatives up to order NQ
1240 C are available from the history array YH, described above under
1241 C optional outputs. In DEWSET, YH is identical to the YCUR array,
1242 C extended to NQ + 1 columns with a column length of NYH and scale
1243 C factors of H**j/factorial(j). On the first call for the problem,
1244 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1245 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
1246 C can be obtained by including in DEWSET the statements:
1247 C DOUBLE PRECISION RLS
1248 C COMMON /DLS001/ RLS(218),ILS(37)
1249 C NQ = ILS(33)
1250 C NST = ILS(34)
1251 C H = RLS(212)
1252 C Thus, for example, the current value of dy/dt can be obtained as
1253 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
1254 C unnecessary when NST = 0).
1255 C
1256 C (b) DVNORM.
1257 C The following is a real function routine which computes the weighted
1258 C root-mean-square norm of a vector v:
1259 C D = DVNORM (N, V, W)
1260 C where:
1261 C N = the length of the vector,
1262 C V = real array of length N containing the vector,
1263 C W = real array of length N containing weights,
1264 C D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).
1265 C DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
1266 C EWT is as set by Subroutine DEWSET.
1267 C
1268 C If the user supplies this function, it should return a non-negative
1269 C value of DVNORM suitable for use in the error control in DLSODIS.
1270 C None of the arguments should be altered by DVNORM.
1271 C For example, a user-supplied DVNORM routine might:
1272 C -substitute a max-norm of (V(i)*w(I)) for the RMS-norm, or
1273 C -ignore some components of V in the norm, with the effect of
1274 C suppressing the error control on those components of y.
1275 C-----------------------------------------------------------------------
1276 C
1277 C***REVISION HISTORY (YYYYMMDD)
1278 C 19820714 DATE WRITTEN
1279 C 19830812 Major update, based on recent LSODI and LSODES revisions:
1280 C Upgraded MDI in ODRV package: operates on M + M-transpose.
1281 C Numerous revisions in use of work arrays;
1282 C use wordlength ratio LENRAT; added IPISP & LRAT to Common;
1283 C added optional outputs IPIAN/IPJAN;
1284 C Added routine CNTNZU; added NZL and NZU to /LSS001/;
1285 C changed ADJLR call logic; added optional outputs NZL & NZU;
1286 C revised counter initializations; revised PREPI stmt. nos.;
1287 C revised difference quotient increment;
1288 C eliminated block /LSI001/, using IERPJ flag;
1289 C revised STODI logic after PJAC return;
1290 C revised tuning of H change and step attempts in STODI;
1291 C corrections to main prologue and comments throughout.
1292 C 19870320 Corrected jump on test of umax in CDRV routine.
1293 C 20010125 Numerous revisions: corrected comments throughout;
1294 C removed TRET from Common; rewrote EWSET with 4 loops;
1295 C fixed t test in INTDY; added Cray directives in STODI;
1296 C in STODI, fixed DELP init. and logic around PJAC call;
1297 C combined routines to save/restore Common;
1298 C passed LEVEL = 0 in error message calls (except run abort).
1299 C 20010425 Major update: convert source lines to upper case;
1300 C added *DECK lines; changed from 1 to * in dummy dimensions;
1301 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
1302 C renamed routines for uniqueness across single/double prec.;
1303 C converted intrinsic names to generic form;
1304 C removed ILLIN and NTREP (data loaded) from Common;
1305 C removed all 'own' variables from Common;
1306 C changed error messages to quoted strings;
1307 C replaced XERRWV/XERRWD with 1993 revised version;
1308 C converted prologues, comments, error messages to mixed case;
1309 C converted arithmetic IF statements to logical IF statements;
1310 C numerous corrections to prologues and internal comments.
1311 C 20010507 Converted single precision source to double precision.
1312 C 20020502 Corrected declarations in descriptions of user routines.
1313 C 20031021 Fixed address offset bugs in Subroutine DPREPI.
1314 C 20031027 Changed 0. to 0.0D0 in Subroutine DPREPI.
1315 C 20031105 Restored 'own' variables to Common blocks, to enable
1316 C interrupt/restart feature.
1317 C 20031112 Added SAVE statements for data-loaded constants.
1318 C 20031117 Changed internal names NRE, LSAVR to NFE, LSAVF resp.
1319 C
1320 C-----------------------------------------------------------------------
1321 C Other routines in the DLSODIS package.
1322 C
1323 C In addition to Subroutine DLSODIS, the DLSODIS package includes the
1324 C following subroutines and function routines:
1325 C DIPREPI acts as an interface between DLSODIS and DPREPI, and also
1326 C does adjusting of work space pointers and work arrays.
1327 C DPREPI is called by DIPREPI to compute sparsity and do sparse
1328 C matrix preprocessing.
1329 C DAINVGS computes the initial value of the vector
1330 C dy/dt = A-inverse * g
1331 C ADJLR adjusts the length of required sparse matrix work space.
1332 C It is called by DPREPI.
1333 C CNTNZU is called by DPREPI and counts the nonzero elements in the
1334 C strict upper triangle of P + P-transpose.
1335 C JGROUP is called by DPREPI to compute groups of Jacobian column
1336 C indices for use when MITER = 2.
1337 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1338 C DSTODI is the core integrator, which does one step of the
1339 C integration and the associated error control.
1340 C DCFODE sets all method coefficients and test constants.
1341 C DPRJIS computes and preprocesses the Jacobian matrix J = dr/dy
1342 C and the Newton iteration matrix P = A - h*l0*J.
1343 C DSOLSS manages solution of linear system in chord iteration.
1344 C DEWSET sets the error weight vector EWT before each step.
1345 C DVNORM computes the weighted RMS-norm of a vector.
1346 C DSRCMS is a user-callable routine to save and restore
1347 C the contents of the internal Common blocks.
1348 C ODRV constructs a reordering of the rows and columns of
1349 C a matrix by the minimum degree algorithm. ODRV is a
1350 C driver routine which calls Subroutines MD, MDI, MDM,
1351 C MDP, MDU, and SRO. See Ref. 2 for details. (The ODRV
1352 C module has been modified since Ref. 2, however.)
1353 C CDRV performs reordering, symbolic factorization, numerical
1354 C factorization, or linear system solution operations,
1355 C depending on a path argument IPATH. CDRV is a
1356 C driver routine which calls Subroutines NROC, NSFC,
1357 C NNFC, NNSC, and NNTC. See Ref. 3 for details.
1358 C DLSODIS uses CDRV to solve linear systems in which the
1359 C coefficient matrix is P = A - con*J, where A is the
1360 C matrix for the linear system A(t,y)*dy/dt = g(t,y),
1361 C con is a scalar, and J is an approximation to
1362 C the Jacobian dr/dy. Because CDRV deals with rowwise
1363 C sparsity descriptions, CDRV works with P-transpose, not P.
1364 C DLSODIS also uses CDRV to solve the linear system
1365 C A(t,y)*dy/dt = g(t,y) for dy/dt when ISTATE = 0.
1366 C (For this, CDRV works with A-transpose, not A.)
1367 C DUMACH computes the unit roundoff in a machine-independent manner.
1368 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
1369 C error messages and warnings. XERRWD is machine-dependent.
1370 C Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1371 C All the others are subroutines.
1372 C
1373 C-----------------------------------------------------------------------
1394 LOGICAL IHIT
1395 CHARACTER*60 MSG
1397 C-----------------------------------------------------------------------
1398 C The following two internal Common blocks contain
1399 C (a) variables which are local to any subroutine but whose values must
1400 C be preserved between calls to the routine ("own" variables), and
1401 C (b) variables which are communicated between subroutines.
1402 C The block DLS001 is declared in subroutines DLSODIS, DIPREPI, DPREPI,
1403 C DINTDY, DSTODI, DPRJIS, and DSOLSS.
1404 C The block DLSS01 is declared in subroutines DLSODIS, DAINVGS,
1405 C DIPREPI, DPREPI, DPRJIS, and DSOLSS.
1406 C Groups of variables are replaced by dummy arrays in the Common
1407 C declarations in routines where those variables are not used.
1408 C-----------------------------------------------------------------------
1415 C
1421 C
1423 C-----------------------------------------------------------------------
1424 C In the Data statement below, set LENRAT equal to the ratio of
1425 C the wordlength for a real number to that for an integer. Usually,
1426 C LENRAT = 1 for single precision and 2 for double precision. If the
1427 C true ratio is not an integer, use the next smaller integer (.ge. 1),
1428 C-----------------------------------------------------------------------
1430 C-----------------------------------------------------------------------
1431 C Block A.
1432 C This code block is executed on every call.
1433 C It tests ISTATE and ITASK for legality and branches appropirately.
1434 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1435 C not yet been done, an error return occurs.
1436 C If ISTATE = 0 or 1 and TOUT = T, return immediately.
1437 C-----------------------------------------------------------------------
1446 C-----------------------------------------------------------------------
1447 C Block B.
1448 C The next code block is executed for the initial call (ISTATE = 0 or 1)
1449 C or for a continuation call with parameter changes (ISTATE = 3).
1450 C It contains checking of all inputs and various initializations.
1451 C If ISTATE = 0 or 1, the final setting of work space pointers, the
1452 C matrix preprocessing, and other initializations are done in Block C.
1453 C
1454 C First check legality of the non-optional inputs NEQ, ITOL, IOPT, and
1455 C MF.
1456 C-----------------------------------------------------------------------
1474 C Next process and check the optional inputs. --------------------------
1477 MXSTEP = MXSTP0
1478 MXHNIL = MXHNL0
1502 C Check RTOL and ATOL for legality. ------------------------------------
1511 C-----------------------------------------------------------------------
1512 C Compute required work array lengths, as far as possible, and test
1513 C these against LRW and LIW. Then set tentative pointers for work
1514 C arrays. Pointers to RWORK/IWORK segments are named by prefixing L to
1515 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1516 C Segments of RWORK (in order) are denoted WM, YH, SAVR, EWT, ACOR.
1517 C The required length of the matrix work space WM is not yet known,
1518 C and so a crude minimum value is used for the initial tests of LRW
1519 C and LIW, and YH is temporarily stored as far to the right in RWORK
1520 C as possible, to leave the maximum amount of space for WM for matrix
1521 C preprocessing. Thus if MOSS .ne. 2 or 4, some of the segments of
1522 C RWORK are temporarily omitted, as they are not needed in the
1523 C preprocessing. These omitted segments are: ACOR if ISTATE = 1,
1524 C EWT and ACOR if ISTATE = 3 and MOSS = 1, and SAVR, EWT, and ACOR if
1525 C ISTATE = 3 and MOSS = 0.
1526 C-----------------------------------------------------------------------
1527 LRAT = LENRAT
1561 LENYHM = NCOLM*NYH
1562 LENYHT = LENYHM
1569 LENWK = LWTEM
1574 ISTATC = ISTATE
1576 C-----------------------------------------------------------------------
1577 C ISTATE = 3. Move YH to its new location.
1578 C Note that only the part of YH needed for the next step, namely
1579 C MIN(NQ+1,MAXORD+2) columns, is actually moved.
1580 C A temporary error weight array EWT is loaded if MOSS = 2 or 4.
1581 C Sparse matrix processing is done in DIPREPI/DPREPI.
1582 C If MAXORD was reduced below NQ, then the pointers are finally set
1583 C so that SAVR is identical to (YH*,MAXORD+2)
1584 C-----------------------------------------------------------------------
1587 C Move YH. Move right if LYHD < 0; move left if LYHD > 0. -------------
1592 ENDIF
1596 ENDIF
1600 C Temporarily load EWT if MOSS = 2 or 4.
1606 C DIPREPI and DPREPI do sparse matrix preprocessing. -------------------
1621 C Set flag to signal changes to DSTODI.---------------------------------
1624 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into YDOTI. --------
1628 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1635 C-----------------------------------------------------------------------
1636 C Block C.
1637 C The next block is for the initial call only (ISTATE = 0 or 1).
1638 C It contains all remaining initializations, the call to DAINVGS
1639 C (if ISTATE = 0), the sparse matrix preprocessing, and the
1640 C calculation if the initial step size.
1641 C The error weights in EWT are inverted after being loaded.
1642 C-----------------------------------------------------------------------
1644 LYH = LYHN
1646 TN = T
1654 C Load the initial value vector in YH.----------------------------------
1658 C Initial dy/dt was supplied. Load it into YH (LYD0 points to YH(*,2).)
1663 C Load and invert the EWT array. (H is temporarily set to 1.0.)--------
1668 C Call DIPREPI and DPREPI to do sparse matrix preprocessing.------------
1680 C Compute initial dy/dt, if necessary, and load it into YH.-------------
1688 C Check TCRIT for legality (ITASK = 4 or 5). ---------------------------
1695 C Initialize all remaining parameters. ---------------------------------
1709 C-----------------------------------------------------------------------
1710 C The coding below computes the step size, H0, to be attempted on the
1711 C first step, unless the user has supplied a value for this.
1712 C First check that TOUT - T differs significantly from zero.
1713 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1714 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1715 C so as to be between 100*UROUND and 1.0E-3.
1716 C Then the computed value H0 is given by..
1717 C NEQ
1718 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * Sum ( YDOT(i)/ywt(i) )**2 )
1719 C 1
1720 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1721 C YDOT(i) = i-th component of initial value of dy/dt,
1722 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1723 C The sign of H0 is inferred from the initial values of TOUT and T.
1724 C-----------------------------------------------------------------------
1747 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1750 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1751 H = H0
1755 C-----------------------------------------------------------------------
1756 C Block D.
1757 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1758 C and is to check stop conditions before taking a step.
1759 C-----------------------------------------------------------------------
1765 T = TOUT
1777 T = TOUT
1788 C-----------------------------------------------------------------------
1789 C Block E.
1790 C The next block is normally executed for all calls and contains
1791 C the call to the one-step core integrator DSTODI.
1792 C
1793 C This is a looping point for the integration steps.
1794 C
1795 C First check for too many steps being taken, update EWT (if not at
1796 C start of problem), check for too much accuracy being requested, and
1797 C check for H below the roundoff level in T.
1798 C-----------------------------------------------------------------------
1825 C-----------------------------------------------------------------------
1826 C CALL DSTODI(NEQ,Y,YH,NYH,YH1,EWT,SAVF,SAVR,ACOR,WM,WM,RES,
1827 C ADDA,JAC,DPRJIS,DSOLSS)
1828 C Note: SAVF in DSTODI occupies the same space as YDOTI in DLSODIS.
1829 C-----------------------------------------------------------------------
1835 C
1836 C KGO = 1:success; 2:error test failure; 3:convergence failure;
1837 C 4:RES ordered return; 5:RES returned error;
1838 C 6:fatal error from CDRV via DPRJIS or DSOLSS.
1839 C-----------------------------------------------------------------------
1840 C Block F.
1841 C The following block handles the case of a successful return from the
1842 C core integrator (KFLAG = 0). Test for stop conditions.
1843 C-----------------------------------------------------------------------
1846 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1849 T = TOUT
1851 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1854 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1857 T = TOUT
1867 C ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1870 C-----------------------------------------------------------------------
1871 C Block G.
1872 C The following block handles all successful returns from DLSODIS.
1873 C if ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1874 C ISTATE is set to 2, and the optional outputs are loaded into the
1875 C work arrays before returning.
1876 C-----------------------------------------------------------------------
1879 T = TN
1897 RETURN
1898 C-----------------------------------------------------------------------
1899 C Block H.
1900 C The following block handles all unsuccessful returns other than
1901 C those for illegal input. First the error message routine is called.
1902 C If there was an error test or convergence test failure, IMXER is set.
1903 C Then Y is loaded from YH and T is set to TN.
1904 C The optional outputs are loaded into the work arrays before returning.
1905 C-----------------------------------------------------------------------
1906 C The maximum number of steps was taken before reaching TOUT. ----------
1913 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1919 C Too much accuracy requested for machine precision. -------------------
1927 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1934 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1943 C IRES = 3 returned by RES, despite retries by DSTODI. -----------------
1950 C KFLAG = -5. Fatal error flag returned by DPRJIS or DSOLSS (CDRV). ---
1959 C DAINVGS failed because matrix A was singular. ------------------------
1967 RETURN
1968 C DAINVGS failed because RES set IRES to 2 or 3. -----------------------
1976 RETURN
1977 C Compute IMXER if relevant. -------------------------------------------
1983 BIG = SIZE
1984 IMXER = I
1987 C Compute residual if relevant. ----------------------------------------
2001 C set y vector, t, and optional outputs. -------------------------------
2018 RETURN
2019 C-----------------------------------------------------------------------
2020 C Block I.
2021 C The following block handles all error returns due to illegal input
2022 C (ISTATE = -3), as detected before calling the core integrator.
2023 C First the error message routine is called. If the illegal input
2024 C is a negative ISTATE, the run is aborted (apparent infinite loop).
2025 C-----------------------------------------------------------------------
2148 ENDIF
2152 ENDIF
2159 C
2161 RETURN
2162 C
2165 RETURN
2166 C----------------------- End of Subroutine DLSODIS ---------------------
2167 END