| blob | dbd7629ac7b3931496af96f82a27a344bb7983d2 |
1 *DECK DLSOIBT
9 C-----------------------------------------------------------------------
10 C This is the 18 November 2003 version of
11 C DLSOIBT: Livermore Solver for Ordinary differential equations given
12 C in Implicit form, with Block-Tridiagonal Jacobian treatment.
13 C
14 C This version is in double precision.
15 C
16 C DLSOIBT 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 DLSOIBT is a variant version of the DLSODI package, for the case where
29 C the matrices A, dg/dy, and d(A*s)/dy are all block-tridiagonal.
30 C-----------------------------------------------------------------------
31 C Reference:
32 C Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE
33 C Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.),
34 C North-Holland, Amsterdam, 1983, pp. 55-64.
35 C-----------------------------------------------------------------------
36 C Authors: Alan C. Hindmarsh and Jeffrey F. Painter
37 C Center for Applied Scientific Computing, L-561
38 C Lawrence Livermore National Laboratory
39 C Livermore, CA 94551
40 C and
41 C Charles S. Kenney
42 C formerly at: Naval Weapons Center
43 C China Lake, CA 93555
44 C-----------------------------------------------------------------------
45 C Summary of Usage.
46 C
47 C Communication between the user and the DLSOIBT package, for normal
48 C situations, is summarized here. This summary describes only a subset
49 C of the full set of options available. See the full description for
50 C details, including optional communication, nonstandard options,
51 C and instructions for special situations. See also the example
52 C problem (with program and output) following this summary.
53 C
54 C A. First, provide a subroutine of the form:
55 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
56 C DOUBLE PRECISION T, Y(*), S(*), R(*)
57 C which computes the residual function
58 C r = g(t,y) - A(t,y) * s ,
59 C as a function of t and the vectors y and s. (s is an internally
60 C generated approximation to dy/dt.) The arrays Y and S are inputs
61 C to the RES routine and should not be altered. The residual
62 C vector is to be stored in the array R. The argument IRES should be
63 C ignored for casual use of DLSOIBT. (For uses of IRES, see the
64 C paragraph on RES in the full description below.)
65 C
66 C B. Next, identify the block structure of the matrices A = A(t,y) and
67 C dr/dy. DLSOIBT must deal internally with a linear combination, P, of
68 C these two matrices. The matrix P (hence both A and dr/dy) must have
69 C a block-tridiagonal form with fixed structure parameters
70 C MB = block size, MB .ge. 1, and
71 C NB = number of blocks in each direction, NB .ge. 4,
72 C with MB*NB = NEQ. In each of the NB block-rows of the matrix P
73 C (each consisting of MB consecutive rows), the nonzero elements are
74 C to lie in three consecutive MB by MB blocks. In block-rows
75 C 2 through NB - 1, these are centered about the main diagonal.
76 C in block-rows 1 and NB, they are the diagonal blocks and the two
77 C blocks adjacent to the diagonal block. (Thus block positions (1,3)
78 C and (NB,NB-2) can be nonzero.)
79 C Alternatively, P (hence A and dr/dy) may be only approximately
80 C equal to matrices with this form, and DLSOIBT should still succeed.
81 C The block-tridiagonal matrix P is described by three arrays,
82 C each of size MB by MB by NB:
83 C PA = array of diagonal blocks,
84 C PB = array of superdiagonal (and one subdiagonal) blocks, and
85 C PC = array of subdiagonal (and one superdiagonal) blocks.
86 C Specifically, the three MB by MB blocks in the k-th block-row of P
87 C are stored in (reading across):
88 C PC(*,*,k) = block to the left of the diagonal block,
89 C PA(*,*,k) = diagonal block, and
90 C PB(*,*,k) = block to the right of the diagonal block,
91 C except for k = 1, where the three blocks (reading across) are
92 C PA(*,*,1) (= diagonal block), PB(*,*,1), and PC(*,*,1),
93 C and k = NB, where they are
94 C PB(*,*,NB), PC(*,*,NB), and PA(*,*,NB) (= diagonal block).
95 C (Each asterisk * stands for an index that ranges from 1 to MB.)
96 C
97 C C. You must also provide a subroutine of the form:
98 C SUBROUTINE ADDA (NEQ, T, Y, MB, NB, PA, PB, PC)
99 C DOUBLE PRECISION T, Y(*), PA(MB,MB,NB), PB(MB,MB,NB), PC(MB,MB,NB)
100 C which adds the nonzero blocks of the matrix A = A(t,y) to the
101 C contents of the arrays PA, PB, and PC, following the structure
102 C description in Paragraph B above.
103 C T and the Y array are input and should not be altered.
104 C Thus the affect of ADDA should be the following:
105 C DO 30 K = 1,NB
106 C DO 20 J = 1,MB
107 C DO 10 I = 1,MB
108 C PA(I,J,K) = PA(I,J,K) +
109 C ( (I,J) element of K-th diagonal block of A)
110 C PB(I,J,K) = PB(I,J,K) +
111 C ( (I,J) element of block in block position (K,K+1) of A,
112 C or in block position (NB,NB-2) if K = NB)
113 C PC(I,J,K) = PC(I,J,K) +
114 C ( (I,J) element of block in block position (K,K-1) of A,
115 C or in block position (1,3) if K = 1)
116 C 10 CONTINUE
117 C 20 CONTINUE
118 C 30 CONTINUE
119 C
120 C D. For the sake of efficiency, you are encouraged to supply the
121 C Jacobian matrix dr/dy in closed form, where r = g(t,y) - A(t,y)*s
122 C (s = a fixed vector) as above. If dr/dy is being supplied,
123 C use MF = 21, and provide a subroutine of the form:
124 C SUBROUTINE JAC (NEQ, T, Y, S, MB, NB, PA, PB, PC)
125 C DOUBLE PRECISION T, Y(*), S(*), PA(MB,MB,NB), PB(MB,MB,NB),
126 C 1 PC(MB,MB,NB)
127 C which computes dr/dy as a function of t, y, and s. Here T, Y, and
128 C S are inputs, and the routine is to load dr/dy into PA, PB, PC,
129 C according to the structure description in Paragraph B above.
130 C That is, load the diagonal blocks into PA, the superdiagonal blocks
131 C (and block (NB,NB-2) ) into PB, and the subdiagonal blocks (and
132 C block (1,3) ) into PC. The blocks in block-row k of dr/dy are to
133 C be loaded into PA(*,*,k), PB(*,*,k), and PC(*,*,k).
134 C Only nonzero elements need be loaded, and the indexing
135 C of PA, PB, and PC is the same as in the ADDA routine.
136 C Note that if A is independent of Y (or this dependence
137 C is weak enough to be ignored) then JAC is to compute dg/dy.
138 C If it is not feasible to provide a JAC routine, use
139 C MF = 22, and DLSOIBT will compute an approximate Jacobian
140 C internally by difference quotients.
141 C
142 C E. Next decide whether or not to provide the initial value of the
143 C derivative vector dy/dt. If the initial value of A(t,y) is
144 C nonsingular (and not too ill-conditioned), you may let DLSOIBT compute
145 C this vector (ISTATE = 0). (DLSOIBT will solve the system A*s = g for
146 C s, with initial values of A and g.) If A(t,y) is initially
147 C singular, then the system is a differential-algebraic system, and
148 C you must make use of the particular form of the system to compute the
149 C initial values of y and dy/dt. In that case, use ISTATE = 1 and
150 C load the initial value of dy/dt into the array YDOTI.
151 C The input array YDOTI and the initial Y array must be consistent with
152 C the equations A*dy/dt = g. This implies that the initial residual
153 C r = g(t,y) - A(t,y)*YDOTI must be approximately zero.
154 C
155 C F. Write a main program which calls Subroutine DLSOIBT once for
156 C each point at which answers are desired. This should also provide
157 C for possible use of logical unit 6 for output of error messages by
158 C DLSOIBT. on the first call to DLSOIBT, supply arguments as follows:
159 C RES = name of user subroutine for residual function r.
160 C ADDA = name of user subroutine for computing and adding A(t,y).
161 C JAC = name of user subroutine for Jacobian matrix dr/dy
162 C (MF = 21). If not used, pass a dummy name.
163 C Note: the names for the RES and ADDA routines and (if used) the
164 C JAC routine must be declared External in the calling program.
165 C NEQ = number of scalar equations in the system.
166 C Y = array of initial values, of length NEQ.
167 C YDOTI = array of length NEQ (containing initial dy/dt if ISTATE = 1).
168 C T = the initial value of the independent variable.
169 C TOUT = first point where output is desired (.ne. T).
170 C ITOL = 1 or 2 according as ATOL (below) is a scalar or array.
171 C RTOL = relative tolerance parameter (scalar).
172 C ATOL = absolute tolerance parameter (scalar or array).
173 C the estimated local error in y(i) will be controlled so as
174 C to be roughly less (in magnitude) than
175 C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or
176 C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2.
177 C Thus the local error test passes if, in each component,
178 C either the absolute error is less than ATOL (or ATOL(i)),
179 C or the relative error is less than RTOL.
180 C Use RTOL = 0.0 for pure absolute error control, and
181 C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error
182 C control. Caution: Actual (global) errors may exceed these
183 C local tolerances, so choose them conservatively.
184 C ITASK = 1 for normal computation of output values of y at t = TOUT.
185 C ISTATE = integer flag (input and output). Set ISTATE = 1 if the
186 C initial dy/dt is supplied, and 0 otherwise.
187 C IOPT = 0 to indicate no optional inputs used.
188 C RWORK = real work array of length at least:
189 C 22 + 9*NEQ + 3*MB*MB*NB for MF = 21 or 22.
190 C LRW = declared length of RWORK (in user's dimension).
191 C IWORK = integer work array of length at least 20 + NEQ.
192 C Input in IWORK(1) the block size MB and in IWORK(2) the
193 C number NB of blocks in each direction along the matrix A.
194 C These must satisfy MB .ge. 1, NB .ge. 4, and MB*NB = NEQ.
195 C LIW = declared length of IWORK (in user's dimension).
196 C MF = method flag. Standard values are:
197 C 21 for a user-supplied Jacobian.
198 C 22 for an internally generated Jacobian.
199 C For other choices of MF, see the paragraph on MF in
200 C the full description below.
201 C Note that the main program must declare arrays Y, YDOTI, RWORK, IWORK,
202 C and possibly ATOL.
203 C
204 C G. The output from the first call (or any call) is:
205 C Y = array of computed values of y(t) vector.
206 C T = corresponding value of independent variable (normally TOUT).
207 C ISTATE = 2 if DLSOIBT was successful, negative otherwise.
208 C -1 means excess work done on this call (check all inputs).
209 C -2 means excess accuracy requested (tolerances too small).
210 C -3 means illegal input detected (see printed message).
211 C -4 means repeated error test failures (check all inputs).
212 C -5 means repeated convergence failures (perhaps bad Jacobian
213 C supplied or wrong choice of tolerances).
214 C -6 means error weight became zero during problem. (Solution
215 C component i vanished, and ATOL or ATOL(i) = 0.)
216 C -7 cannot occur in casual use.
217 C -8 means DLSOIBT was unable to compute the initial dy/dt.
218 C In casual use, this means A(t,y) is initially singular.
219 C Supply YDOTI and use ISTATE = 1 on the first call.
220 C
221 C If DLSOIBT returns ISTATE = -1, -4, or -5, then the output of
222 C DLSOIBT also includes YDOTI = array containing residual vector
223 C r = g - A * dy/dt evaluated at the current t, y, and dy/dt.
224 C
225 C H. To continue the integration after a successful return, simply
226 C reset TOUT and call DLSOIBT again. No other parameters need be reset.
227 C
228 C-----------------------------------------------------------------------
229 C Example Problem.
230 C
231 C The following is an example problem, with the coding needed
232 C for its solution by DLSOIBT. The problem comes from the partial
233 C differential equation (the Burgers equation)
234 C du/dt = - u * du/dx + eta * d**2 u/dx**2, eta = .05,
235 C on -1 .le. x .le. 1. The boundary conditions are
236 C du/dx = 0 at x = -1 and at x = 1.
237 C The initial profile is a square wave,
238 C u = 1 in ABS(x) .lt. .5, u = .5 at ABS(x) = .5, u = 0 elsewhere.
239 C The PDE is discretized in x by a simplified Galerkin method,
240 C using piecewise linear basis functions, on a grid of 40 intervals.
241 C The equations at x = -1 and 1 use a 3-point difference approximation
242 C for the right-hand side. The result is a system A * dy/dt = g(y),
243 C of size NEQ = 41, where y(i) is the approximation to u at x = x(i),
244 C with x(i) = -1 + (i-1)*delx, delx = 2/(NEQ-1) = .05. The individual
245 C equations in the system are
246 C dy(1)/dt = ( y(3) - 2*y(2) + y(1) ) * eta / delx**2,
247 C dy(NEQ)/dt = ( y(NEQ-2) - 2*y(NEQ-1) + y(NEQ) ) * eta / delx**2,
248 C and for i = 2, 3, ..., NEQ-1,
249 C (1/6) dy(i-1)/dt + (4/6) dy(i)/dt + (1/6) dy(i+1)/dt
250 C = ( y(i-1)**2 - y(i+1)**2 ) / (4*delx)
251 C + ( y(i+1) - 2*y(i) + y(i-1) ) * eta / delx**2.
252 C The following coding solves the problem with MF = 21, with output
253 C of solution statistics at t = .1, .2, .3, and .4, and of the
254 C solution vector at t = .4. Here the block size is just MB = 1.
255 C
256 C EXTERNAL RESID, ADDABT, JACBT
257 C DOUBLE PRECISION ATOL, RTOL, RWORK, T, TOUT, Y, YDOTI
258 C DIMENSION Y(41), YDOTI(41), RWORK(514), IWORK(61)
259 C NEQ = 41
260 C DO 10 I = 1,NEQ
261 C 10 Y(I) = 0.0
262 C Y(11) = 0.5
263 C DO 20 I = 12,30
264 C 20 Y(I) = 1.0
265 C Y(31) = 0.5
266 C T = 0.0
267 C TOUT = 0.1
268 C ITOL = 1
269 C RTOL = 1.0D-4
270 C ATOL = 1.0D-5
271 C ITASK = 1
272 C ISTATE = 0
273 C IOPT = 0
274 C LRW = 514
275 C LIW = 61
276 C IWORK(1) = 1
277 C IWORK(2) = NEQ
278 C MF = 21
279 C DO 40 IO = 1,4
280 C CALL DLSOIBT (RESID, ADDABT, JACBT, NEQ, Y, YDOTI, T, TOUT,
281 C 1 ITOL,RTOL,ATOL, ITASK, ISTATE, IOPT, RWORK,LRW,IWORK,LIW, MF)
282 C WRITE (6,30) T, IWORK(11), IWORK(12), IWORK(13)
283 C 30 FORMAT(' At t =',F5.2,' No. steps =',I4,' No. r-s =',I4,
284 C 1 ' No. J-s =',I3)
285 C IF (ISTATE .NE. 2) GO TO 90
286 C TOUT = TOUT + 0.1
287 C 40 CONTINUE
288 C WRITE(6,50) (Y(I),I=1,NEQ)
289 C 50 FORMAT(/' Final solution values..'/9(5D12.4/))
290 C STOP
291 C 90 WRITE(6,95) ISTATE
292 C 95 FORMAT(///' Error halt.. ISTATE =',I3)
293 C STOP
294 C END
295 C
296 C SUBROUTINE RESID (N, T, Y, S, R, IRES)
297 C DOUBLE PRECISION T, Y, S, R, ETA, DELX, EODSQ
298 C DIMENSION Y(N), S(N), R(N)
299 C DATA ETA/0.05/, DELX/0.05/
300 C EODSQ = ETA/DELX**2
301 C R(1) = EODSQ*(Y(3) - 2.0*Y(2) + Y(1)) - S(1)
302 C NM1 = N - 1
303 C DO 10 I = 2,NM1
304 C R(I) = (Y(I-1)**2 - Y(I+1)**2)/(4.0*DELX)
305 C 1 + EODSQ*(Y(I+1) - 2.0*Y(I) + Y(I-1))
306 C 2 - (S(I-1) + 4.0*S(I) + S(I+1))/6.0
307 C 10 CONTINUE
308 C R(N) = EODSQ*(Y(N-2) - 2.0*Y(NM1) + Y(N)) - S(N)
309 C RETURN
310 C END
311 C
312 C SUBROUTINE ADDABT (N, T, Y, MB, NB, PA, PB, PC)
313 C DOUBLE PRECISION T, Y, PA, PB, PC
314 C DIMENSION Y(N), PA(MB,MB,NB), PB(MB,MB,NB), PC(MB,MB,NB)
315 C PA(1,1,1) = PA(1,1,1) + 1.0
316 C NM1 = N - 1
317 C DO 10 K = 2,NM1
318 C PA(1,1,K) = PA(1,1,K) + (4.0/6.0)
319 C PB(1,1,K) = PB(1,1,K) + (1.0/6.0)
320 C PC(1,1,K) = PC(1,1,K) + (1.0/6.0)
321 C 10 CONTINUE
322 C PA(1,1,N) = PA(1,1,N) + 1.0
323 C RETURN
324 C END
325 C
326 C SUBROUTINE JACBT (N, T, Y, S, MB, NB, PA, PB, PC)
327 C DOUBLE PRECISION T, Y, S, PA, PB, PC, ETA, DELX, EODSQ
328 C DIMENSION Y(N), S(N), PA(MB,MB,NB),PB(MB,MB,NB),PC(MB,MB,NB)
329 C DATA ETA/0.05/, DELX/0.05/
330 C EODSQ = ETA/DELX**2
331 C PA(1,1,1) = EODSQ
332 C PB(1,1,1) = -2.0*EODSQ
333 C PC(1,1,1) = EODSQ
334 C DO 10 K = 2,N
335 C PA(1,1,K) = -2.0*EODSQ
336 C PB(1,1,K) = -Y(K+1)*(0.5/DELX) + EODSQ
337 C PC(1,1,K) = Y(K-1)*(0.5/DELX) + EODSQ
338 C 10 CONTINUE
339 C PB(1,1,N) = EODSQ
340 C PC(1,1,N) = -2.0*EODSQ
341 C PA(1,1,N) = EODSQ
342 C RETURN
343 C END
344 C
345 C The output of this program (on a CDC-7600 in single precision)
346 C is as follows:
347 C
348 C At t = 0.10 No. steps = 35 No. r-s = 45 No. J-s = 9
349 C At t = 0.20 No. steps = 43 No. r-s = 54 No. J-s = 10
350 C At t = 0.30 No. steps = 48 No. r-s = 60 No. J-s = 11
351 C At t = 0.40 No. steps = 51 No. r-s = 64 No. J-s = 12
352 C
353 C Final solution values..
354 C 1.2747e-02 1.1997e-02 1.5560e-02 2.3767e-02 3.7224e-02
355 C 5.6646e-02 8.2645e-02 1.1557e-01 1.5541e-01 2.0177e-01
356 C 2.5397e-01 3.1104e-01 3.7189e-01 4.3530e-01 5.0000e-01
357 C 5.6472e-01 6.2816e-01 6.8903e-01 7.4612e-01 7.9829e-01
358 C 8.4460e-01 8.8438e-01 9.1727e-01 9.4330e-01 9.6281e-01
359 C 9.7632e-01 9.8426e-01 9.8648e-01 9.8162e-01 9.6617e-01
360 C 9.3374e-01 8.7535e-01 7.8236e-01 6.5321e-01 5.0003e-01
361 C 3.4709e-01 2.1876e-01 1.2771e-01 7.3671e-02 5.0642e-02
362 C 5.4496e-02
363 C
364 C-----------------------------------------------------------------------
365 C Full Description of User Interface to DLSOIBT.
366 C
367 C The user interface to DLSOIBT consists of the following parts.
368 C
369 C 1. The call sequence to Subroutine DLSOIBT, which is a driver
370 C routine for the solver. This includes descriptions of both
371 C the call sequence arguments and of user-supplied routines.
372 C Following these descriptions is a description of
373 C optional inputs available through the call sequence, and then
374 C a description of optional outputs (in the work arrays).
375 C
376 C 2. Descriptions of other routines in the DLSOIBT package that may be
377 C (optionally) called by the user. These provide the ability to
378 C alter error message handling, save and restore the internal
379 C Common, and obtain specified derivatives of the solution y(t).
380 C
381 C 3. Descriptions of Common blocks to be declared in overlay
382 C or similar environments, or to be saved when doing an interrupt
383 C of the problem and continued solution later.
384 C
385 C 4. Description of two routines in the DLSOIBT package, either of
386 C which the user may replace with his/her own version, if desired.
387 C These relate to the measurement of errors.
388 C
389 C-----------------------------------------------------------------------
390 C Part 1. Call Sequence.
391 C
392 C The call sequence parameters used for input only are
393 C RES, ADDA, JAC, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK,
394 C IOPT, LRW, LIW, MF,
395 C and those used for both input and output are
396 C Y, T, ISTATE, YDOTI.
397 C The work arrays RWORK and IWORK are also used for additional and
398 C optional inputs and optional outputs. (The term output here refers
399 C to the return from Subroutine DLSOIBT to the user's calling program.)
400 C
401 C The legality of input parameters will be thoroughly checked on the
402 C initial call for the problem, but not checked thereafter unless a
403 C change in input parameters is flagged by ISTATE = 3 on input.
404 C
405 C The descriptions of the call arguments are as follows.
406 C
407 C RES = the name of the user-supplied subroutine which supplies
408 C the residual vector for the ODE system, defined by
409 C r = g(t,y) - A(t,y) * s
410 C as a function of the scalar t and the vectors
411 C s and y (s approximates dy/dt). This subroutine
412 C is to have the form
413 C SUBROUTINE RES (NEQ, T, Y, S, R, IRES)
414 C DOUBLE PRECISION T, Y(*), S(*), R(*)
415 C where NEQ, T, Y, S, and IRES are input, and R and
416 C IRES are output. Y, S, and R are arrays of length NEQ.
417 C On input, IRES indicates how DLSOIBT will use the
418 C returned array R, as follows:
419 C IRES = 1 means that DLSOIBT needs the full residual,
420 C r = g - A*s, exactly.
421 C IRES = -1 means that DLSOIBT is using R only to compute
422 C the Jacobian dr/dy by difference quotients.
423 C The RES routine can ignore IRES, or it can omit some terms
424 C if IRES = -1. If A does not depend on y, then RES can
425 C just return R = g when IRES = -1. If g - A*s contains other
426 C additive terms that are independent of y, these can also be
427 C dropped, if done consistently, when IRES = -1.
428 C The subroutine should set the flag IRES if it
429 C encounters a halt condition or illegal input.
430 C Otherwise, it should not reset IRES. On output,
431 C IRES = 1 or -1 represents a normal return, and
432 C DLSOIBT continues integrating the ODE. Leave IRES
433 C unchanged from its input value.
434 C IRES = 2 tells DLSOIBT to immediately return control
435 C to the calling program, with ISTATE = 3. This lets
436 C the calling program change parameters of the problem
437 C if necessary.
438 C IRES = 3 represents an error condition (for example, an
439 C illegal value of y). DLSOIBT tries to integrate the system
440 C without getting IRES = 3 from RES. If it cannot, DLSOIBT
441 C returns with ISTATE = -7 or -1.
442 C On an DLSOIBT return with ISTATE = 3, -1, or -7, the
443 C values of T and Y returned correspond to the last point
444 C reached successfully without getting the flag IRES = 2 or 3.
445 C The flag values IRES = 2 and 3 should not be used to
446 C handle switches or root-stop conditions. This is better
447 C done by calling DLSOIBT in a one-step mode and checking the
448 C stopping function for a sign change at each step.
449 C If quantities computed in the RES routine are needed
450 C externally to DLSOIBT, an extra call to RES should be made
451 C for this purpose, for consistent and accurate results.
452 C To get the current dy/dt for the S argument, use DINTDY.
453 C RES must be declared External in the calling
454 C program. See note below for more about RES.
455 C
456 C ADDA = the name of the user-supplied subroutine which adds the
457 C matrix A = A(t,y) to another matrix, P, stored in
458 C block-tridiagonal form. This routine is to have the form
459 C SUBROUTINE ADDA (NEQ, T, Y, MB, NB, PA, PB, PC)
460 C DOUBLE PRECISION T, Y(*), PA(MB,MB,NB), PB(MB,MB,NB),
461 C 1 PC(MB,MB,NB)
462 C where NEQ, T, Y, MB, NB, and the arrays PA, PB, and PC
463 C are input, and the arrays PA, PB, and PC are output.
464 C Y is an array of length NEQ, and the arrays PA, PB, PC
465 C are all MB by MB by NB.
466 C Here a block-tridiagonal structure is assumed for A(t,y),
467 C and also for the matrix P to which A is added here,
468 C as described in Paragraph B of the Summary of Usage above.
469 C Thus the affect of ADDA should be the following:
470 C DO 30 K = 1,NB
471 C DO 20 J = 1,MB
472 C DO 10 I = 1,MB
473 C PA(I,J,K) = PA(I,J,K) +
474 C ( (I,J) element of K-th diagonal block of A)
475 C PB(I,J,K) = PB(I,J,K) +
476 C ( (I,J) element of block (K,K+1) of A,
477 C or block (NB,NB-2) if K = NB)
478 C PC(I,J,K) = PC(I,J,K) +
479 C ( (I,J) element of block (K,K-1) of A,
480 C or block (1,3) if K = 1)
481 C 10 CONTINUE
482 C 20 CONTINUE
483 C 30 CONTINUE
484 C ADDA must be declared External in the calling program.
485 C See note below for more information about ADDA.
486 C
487 C JAC = the name of the user-supplied subroutine which supplies
488 C the Jacobian matrix, dr/dy, where r = g - A*s. JAC is
489 C required if MITER = 1. Otherwise a dummy name can be
490 C passed. This subroutine is to have the form
491 C SUBROUTINE JAC (NEQ, T, Y, S, MB, NB, PA, PB, PC)
492 C DOUBLE PRECISION T, Y(*), S(*), PA(MB,MB,NB),
493 C 1 PB(MB,MB,NB), PC(MB,MB,NB)
494 C where NEQ, T, Y, S, MB, NB, and the arrays PA, PB, and PC
495 C are input, and the arrays PA, PB, and PC are output.
496 C Y and S are arrays of length NEQ, and the arrays PA, PB, PC
497 C are all MB by MB by NB.
498 C PA, PB, and PC are to be loaded with partial derivatives
499 C (elements of the Jacobian matrix) on output, in terms of the
500 C block-tridiagonal structure assumed, as described
501 C in Paragraph B of the Summary of Usage above.
502 C That is, load the diagonal blocks into PA, the
503 C superdiagonal blocks (and block (NB,NB-2) ) into PB, and
504 C the subdiagonal blocks (and block (1,3) ) into PC.
505 C The blocks in block-row k of dr/dy are to be loaded into
506 C PA(*,*,k), PB(*,*,k), and PC(*,*,k).
507 C Thus the affect of JAC should be the following:
508 C DO 30 K = 1,NB
509 C DO 20 J = 1,MB
510 C DO 10 I = 1,MB
511 C PA(I,J,K) = ( (I,J) element of
512 C K-th diagonal block of dr/dy)
513 C PB(I,J,K) = ( (I,J) element of block (K,K+1)
514 C of dr/dy, or block (NB,NB-2) if K = NB)
515 C PC(I,J,K) = ( (I,J) element of block (K,K-1)
516 C of dr/dy, or block (1,3) if K = 1)
517 C 10 CONTINUE
518 C 20 CONTINUE
519 C 30 CONTINUE
520 C PA, PB, and PC are preset to zero by the solver,
521 C so that only the nonzero elements need be loaded by JAC.
522 C Each call to JAC is preceded by a call to RES with the same
523 C arguments NEQ, T, Y, and S. Thus to gain some efficiency,
524 C intermediate quantities shared by both calculations may be
525 C saved in a user Common block by RES and not recomputed by JAC
526 C if desired. Also, JAC may alter the Y array, if desired.
527 C JAC need not provide dr/dy exactly. A crude
528 C approximation will do, so that DLSOIBT may be used when
529 C A and dr/dy are not really block-tridiagonal, but are close
530 C to matrices that are.
531 C JAC must be declared External in the calling program.
532 C See note below for more about JAC.
533 C
534 C Note on RES, ADDA, and JAC:
535 C These subroutines may access user-defined quantities in
536 C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
537 C (dimensioned in the subroutines) and/or Y has length
538 C exceeding NEQ(1). However, these routines should not alter
539 C NEQ(1), Y(1),...,Y(NEQ) or any other input variables.
540 C See the descriptions of NEQ and Y below.
541 C
542 C NEQ = the size of the system (number of first order ordinary
543 C differential equations or scalar algebraic equations).
544 C Used only for input.
545 C NEQ may be decreased, but not increased, during the problem.
546 C If NEQ is decreased (with ISTATE = 3 on input), the
547 C remaining components of Y should be left undisturbed, if
548 C these are to be accessed in RES, ADDA, or JAC.
549 C
550 C Normally, NEQ is a scalar, and it is generally referred to
551 C as a scalar in this user interface description. However,
552 C NEQ may be an array, with NEQ(1) set to the system size.
553 C (The DLSOIBT package accesses only NEQ(1).) In either case,
554 C this parameter is passed as the NEQ argument in all calls
555 C to RES, ADDA, and JAC. Hence, if it is an array,
556 C locations NEQ(2),... may be used to store other integer data
557 C and pass it to RES, ADDA, or JAC. Each such subroutine
558 C must include NEQ in a Dimension statement in that case.
559 C
560 C Y = a real array for the vector of dependent variables, of
561 C length NEQ or more. Used for both input and output on the
562 C first call (ISTATE = 0 or 1), and only for output on other
563 C calls. On the first call, Y must contain the vector of
564 C initial values. On output, Y contains the computed solution
565 C vector, evaluated at t. If desired, the Y array may be used
566 C for other purposes between calls to the solver.
567 C
568 C This array is passed as the Y argument in all calls to RES,
569 C ADDA, and JAC. Hence its length may exceed NEQ,
570 C and locations Y(NEQ+1),... may be used to store other real
571 C data and pass it to RES, ADDA, or JAC. (The DLSOIBT
572 C package accesses only Y(1),...,Y(NEQ). )
573 C
574 C YDOTI = a real array for the initial value of the vector
575 C dy/dt and for work space, of dimension at least NEQ.
576 C
577 C On input:
578 C If ISTATE = 0 then DLSOIBT will compute the initial value
579 C of dy/dt, if A is nonsingular. Thus YDOTI will
580 C serve only as work space and may have any value.
581 C If ISTATE = 1 then YDOTI must contain the initial value
582 C of dy/dt.
583 C If ISTATE = 2 or 3 (continuation calls) then YDOTI
584 C may have any value.
585 C Note: If the initial value of A is singular, then
586 C DLSOIBT cannot compute the initial value of dy/dt, so
587 C it must be provided in YDOTI, with ISTATE = 1.
588 C
589 C On output, when DLSOIBT terminates abnormally with ISTATE =
590 C -1, -4, or -5, YDOTI will contain the residual
591 C r = g(t,y) - A(t,y)*(dy/dt). If r is large, t is near
592 C its initial value, and YDOTI is supplied with ISTATE = 1,
593 C there may have been an incorrect input value of
594 C YDOTI = dy/dt, or the problem (as given to DLSOIBT)
595 C may not have a solution.
596 C
597 C If desired, the YDOTI array may be used for other
598 C purposes between calls to the solver.
599 C
600 C T = the independent variable. On input, T is used only on the
601 C first call, as the initial point of the integration.
602 C On output, after each call, T is the value at which a
603 C computed solution y is evaluated (usually the same as TOUT).
604 C On an error return, T is the farthest point reached.
605 C
606 C TOUT = the next value of t at which a computed solution is desired.
607 C Used only for input.
608 C
609 C When starting the problem (ISTATE = 0 or 1), TOUT may be
610 C equal to T for one call, then should .ne. T for the next
611 C call. For the initial T, an input value of TOUT .ne. T is
612 C used in order to determine the direction of the integration
613 C (i.e. the algebraic sign of the step sizes) and the rough
614 C scale of the problem. Integration in either direction
615 C (forward or backward in t) is permitted.
616 C
617 C If ITASK = 2 or 5 (one-step modes), TOUT is ignored after
618 C the first call (i.e. the first call with TOUT .ne. T).
619 C Otherwise, TOUT is required on every call.
620 C
621 C If ITASK = 1, 3, or 4, the values of TOUT need not be
622 C monotone, but a value of TOUT which backs up is limited
623 C to the current internal T interval, whose endpoints are
624 C TCUR - HU and TCUR (see optional outputs, below, for
625 C TCUR and HU).
626 C
627 C ITOL = an indicator for the type of error control. See
628 C description below under ATOL. Used only for input.
629 C
630 C RTOL = a relative error tolerance parameter, either a scalar or
631 C an array of length NEQ. See description below under ATOL.
632 C Input only.
633 C
634 C ATOL = an absolute error tolerance parameter, either a scalar or
635 C an array of length NEQ. Input only.
636 C
637 C The input parameters ITOL, RTOL, and ATOL determine
638 C the error control performed by the solver. The solver will
639 C control the vector E = (E(i)) of estimated local errors
640 C in y, according to an inequality of the form
641 C RMS-norm of ( E(i)/EWT(i) ) .le. 1,
642 C where EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
643 C and the RMS-norm (root-mean-square norm) here is
644 C RMS-norm(v) = SQRT(sum v(i)**2 / NEQ). Here EWT = (EWT(i))
645 C is a vector of weights which must always be positive, and
646 C the values of RTOL and ATOL should all be non-negative.
647 C The following table gives the types (scalar/array) of
648 C RTOL and ATOL, and the corresponding form of EWT(i).
649 C
650 C ITOL RTOL ATOL EWT(i)
651 C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL
652 C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i)
653 C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL
654 C 4 array scalar RTOL(i)*ABS(Y(i)) + ATOL(i)
655 C
656 C When either of these parameters is a scalar, it need not
657 C be dimensioned in the user's calling program.
658 C
659 C If none of the above choices (with ITOL, RTOL, and ATOL
660 C fixed throughout the problem) is suitable, more general
661 C error controls can be obtained by substituting
662 C user-supplied routines for the setting of EWT and/or for
663 C the norm calculation. See Part 4 below.
664 C
665 C If global errors are to be estimated by making a repeated
666 C run on the same problem with smaller tolerances, then all
667 C components of RTOL and ATOL (i.e. of EWT) should be scaled
668 C down uniformly.
669 C
670 C ITASK = an index specifying the task to be performed.
671 C Input only. ITASK has the following values and meanings.
672 C 1 means normal computation of output values of y(t) at
673 C t = TOUT (by overshooting and interpolating).
674 C 2 means take one step only and return.
675 C 3 means stop at the first internal mesh point at or
676 C beyond t = TOUT and return.
677 C 4 means normal computation of output values of y(t) at
678 C t = TOUT but without overshooting t = TCRIT.
679 C TCRIT must be input as RWORK(1). TCRIT may be equal to
680 C or beyond TOUT, but not behind it in the direction of
681 C integration. This option is useful if the problem
682 C has a singularity at or beyond t = TCRIT.
683 C 5 means take one step, without passing TCRIT, and return.
684 C TCRIT must be input as RWORK(1).
685 C
686 C Note: If ITASK = 4 or 5 and the solver reaches TCRIT
687 C (within roundoff), it will return T = TCRIT (exactly) to
688 C indicate this (unless ITASK = 4 and TOUT comes before TCRIT,
689 C in which case answers at t = TOUT are returned first).
690 C
691 C ISTATE = an index used for input and output to specify the
692 C state of the calculation.
693 C
694 C On input, the values of ISTATE are as follows.
695 C 0 means this is the first call for the problem, and
696 C DLSOIBT is to compute the initial value of dy/dt
697 C (while doing other initializations). See note below.
698 C 1 means this is the first call for the problem, and
699 C the initial value of dy/dt has been supplied in
700 C YDOTI (DLSOIBT will do other initializations).
701 C See note below.
702 C 2 means this is not the first call, and the calculation
703 C is to continue normally, with no change in any input
704 C parameters except possibly TOUT and ITASK.
705 C (If ITOL, RTOL, and/or ATOL are changed between calls
706 C with ISTATE = 2, the new values will be used but not
707 C tested for legality.)
708 C 3 means this is not the first call, and the
709 C calculation is to continue normally, but with
710 C a change in input parameters other than
711 C TOUT and ITASK. Changes are allowed in
712 C NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, MB, NB,
713 C and any of the optional inputs except H0.
714 C (See IWORK description for MB and NB.)
715 C Note: A preliminary call with TOUT = T is not counted
716 C as a first call here, as no initialization or checking of
717 C input is done. (Such a call is sometimes useful for the
718 C purpose of outputting the initial conditions.)
719 C Thus the first call for which TOUT .ne. T requires
720 C ISTATE = 0 or 1 on input.
721 C
722 C On output, ISTATE has the following values and meanings.
723 C 0 or 1 means nothing was done; TOUT = t and
724 C ISTATE = 0 or 1 on input.
725 C 2 means that the integration was performed successfully.
726 C 3 means that the user-supplied Subroutine RES signalled
727 C DLSOIBT to halt the integration and return (IRES = 2).
728 C Integration as far as T was achieved with no occurrence
729 C of IRES = 2, but this flag was set on attempting the
730 C next step.
731 C -1 means an excessive amount of work (more than MXSTEP
732 C steps) was done on this call, before completing the
733 C requested task, but the integration was otherwise
734 C successful as far as T. (MXSTEP is an optional input
735 C and is normally 500.) To continue, the user may
736 C simply reset ISTATE to a value .gt. 1 and call again
737 C (the excess work step counter will be reset to 0).
738 C In addition, the user may increase MXSTEP to avoid
739 C this error return (see below on optional inputs).
740 C -2 means too much accuracy was requested for the precision
741 C of the machine being used. This was detected before
742 C completing the requested task, but the integration
743 C was successful as far as T. To continue, the tolerance
744 C parameters must be reset, and ISTATE must be set
745 C to 3. The optional output TOLSF may be used for this
746 C purpose. (Note: If this condition is detected before
747 C taking any steps, then an illegal input return
748 C (ISTATE = -3) occurs instead.)
749 C -3 means illegal input was detected, before taking any
750 C integration steps. See written message for details.
751 C Note: If the solver detects an infinite loop of calls
752 C to the solver with illegal input, it will cause
753 C the run to stop.
754 C -4 means there were repeated error test failures on
755 C one attempted step, before completing the requested
756 C task, but the integration was successful as far as T.
757 C The problem may have a singularity, or the input
758 C may be inappropriate.
759 C -5 means there were repeated convergence test failures on
760 C one attempted step, before completing the requested
761 C task, but the integration was successful as far as T.
762 C This may be caused by an inaccurate Jacobian matrix.
763 C -6 means EWT(i) became zero for some i during the
764 C integration. Pure relative error control (ATOL(i) = 0.0)
765 C was requested on a variable which has now vanished.
766 C The integration was successful as far as T.
767 C -7 means that the user-supplied Subroutine RES set
768 C its error flag (IRES = 3) despite repeated tries by
769 C DLSOIBT to avoid that condition.
770 C -8 means that ISTATE was 0 on input but DLSOIBT was unable
771 C to compute the initial value of dy/dt. See the
772 C printed message for details.
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 DLSOIBT 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 real working array (double precision).
792 C The length of RWORK must be at least
793 C 20 + NYH*(MAXORD + 1) + 3*NEQ + LENWM where
794 C NYH = the initial value of NEQ,
795 C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
796 C smaller value is given as an optional input),
797 C LENWM = 3*MB*MB*NB + 2.
798 C (See MF description for the definition of METH.)
799 C Thus if MAXORD has its default value and NEQ is constant,
800 C this length is
801 C 22 + 16*NEQ + 3*MB*MB*NB for MF = 11 or 12,
802 C 22 + 9*NEQ + 3*MB*MB*NB for MF = 21 or 22.
803 C The first 20 words of RWORK are reserved for conditional
804 C and optional inputs and optional outputs.
805 C
806 C The following word in RWORK is a conditional input:
807 C RWORK(1) = TCRIT = critical value of t which the solver
808 C is not to overshoot. Required if ITASK is
809 C 4 or 5, and ignored otherwise. (See ITASK.)
810 C
811 C LRW = the length of the array RWORK, as declared by the user.
812 C (This will be checked by the solver.)
813 C
814 C IWORK = an integer work array. The length of IWORK must be at least
815 C 20 + NEQ . The first few words of IWORK are used for
816 C additional and optional inputs and optional outputs.
817 C
818 C The following 2 words in IWORK are additional required
819 C inputs to DLSOIBT:
820 C IWORK(1) = MB = block size
821 C IWORK(2) = NB = number of blocks in the main diagonal
822 C These must satisfy MB .ge. 1, NB .ge. 4, and MB*NB = NEQ.
823 C
824 C LIW = the length of the array IWORK, as declared by the user.
825 C (This will be checked by the solver.)
826 C
827 C Note: The work arrays must not be altered between calls to DLSOIBT
828 C for the same problem, except possibly for the additional and
829 C optional inputs, and except for the last 3*NEQ words of RWORK.
830 C The latter space is used for internal scratch space, and so is
831 C available for use by the user outside DLSOIBT between calls, if
832 C desired (but not for use by RES, ADDA, or JAC).
833 C
834 C MF = the method flag. used only for input. The legal values of
835 C MF are 11, 12, 21, and 22.
836 C MF has decimal digits METH and MITER: MF = 10*METH + MITER.
837 C METH indicates the basic linear multistep method:
838 C METH = 1 means the implicit Adams method.
839 C METH = 2 means the method based on Backward
840 C Differentiation Formulas (BDFS).
841 C The BDF method is strongly preferred for stiff
842 C problems, while the Adams method is preferred when the
843 C problem is not stiff. If the matrix A(t,y) is
844 C nonsingular, stiffness here can be taken to mean that of
845 C the explicit ODE system dy/dt = A-inverse * g. If A is
846 C singular, the concept of stiffness is not well defined.
847 C If you do not know whether the problem is stiff, we
848 C recommend using METH = 2. If it is stiff, the advantage
849 C of METH = 2 over METH = 1 will be great, while if it is
850 C not stiff, the advantage of METH = 1 will be slight.
851 C If maximum efficiency is important, some experimentation
852 C with METH may be necessary.
853 C MITER indicates the corrector iteration method:
854 C MITER = 1 means chord iteration with a user-supplied
855 C block-tridiagonal Jacobian.
856 C MITER = 2 means chord iteration with an internally
857 C generated (difference quotient) block-
858 C tridiagonal Jacobian approximation, using
859 C 3*MB+1 extra calls to RES per dr/dy evaluation.
860 C If MITER = 1, the user must supply a Subroutine JAC
861 C (the name is arbitrary) as described above under JAC.
862 C For MITER = 2, a dummy argument can be used.
863 C-----------------------------------------------------------------------
864 C Optional Inputs.
865 C
866 C The following is a list of the optional inputs provided for in the
867 C call sequence. (See also Part 2.) For each such input variable,
868 C this table lists its name as used in this documentation, its
869 C location in the call sequence, its meaning, and the default value.
870 C The use of any of these inputs requires IOPT = 1, and in that
871 C case all of these inputs are examined. A value of zero for any
872 C of these optional inputs will cause the default value to be used.
873 C Thus to use a subset of the optional inputs, simply preload
874 C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
875 C then set those of interest to nonzero values.
876 C
877 C Name Location Meaning and Default Value
878 C
879 C H0 RWORK(5) the step size to be attempted on the first step.
880 C The default value is determined by the solver.
881 C
882 C HMAX RWORK(6) the maximum absolute step size allowed.
883 C The default value is infinite.
884 C
885 C HMIN RWORK(7) the minimum absolute step size allowed.
886 C The default value is 0. (This lower bound is not
887 C enforced on the final step before reaching TCRIT
888 C when ITASK = 4 or 5.)
889 C
890 C MAXORD IWORK(5) the maximum order to be allowed. The default
891 C value is 12 if METH = 1, and 5 if METH = 2.
892 C If MAXORD exceeds the default value, it will
893 C be reduced to the default value.
894 C If MAXORD is changed during the problem, it may
895 C cause the current order to be reduced.
896 C
897 C MXSTEP IWORK(6) maximum number of (internally defined) steps
898 C allowed during one call to the solver.
899 C The default value is 500.
900 C
901 C MXHNIL IWORK(7) maximum number of messages printed (per problem)
902 C warning that T + H = T on a step (H = step size).
903 C This must be positive to result in a non-default
904 C value. The default value is 10.
905 C-----------------------------------------------------------------------
906 C Optional Outputs.
907 C
908 C As optional additional output from DLSOIBT, the variables listed
909 C below are quantities related to the performance of DLSOIBT
910 C which are available to the user. These are communicated by way of
911 C the work arrays, but also have internal mnemonic names as shown.
912 C Except where stated otherwise, all of these outputs are defined
913 C on any successful return from DLSOIBT, and on any return with
914 C ISTATE = -1, -2, -4, -5, -6, or -7. On a return with -3 (illegal
915 C input) or -8, they will be unchanged from their existing values
916 C (if any), except possibly for TOLSF, LENRW, and LENIW.
917 C On any error return, outputs relevant to the error will be defined,
918 C as noted below.
919 C
920 C Name Location Meaning
921 C
922 C HU RWORK(11) the step size in t last used (successfully).
923 C
924 C HCUR RWORK(12) the step size to be attempted on the next step.
925 C
926 C TCUR RWORK(13) the current value of the independent variable
927 C which the solver has actually reached, i.e. the
928 C current internal mesh point in t. On output, TCUR
929 C will always be at least as far as the argument
930 C T, but may be farther (if interpolation was done).
931 C
932 C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
933 C computed when a request for too much accuracy was
934 C detected (ISTATE = -3 if detected at the start of
935 C the problem, ISTATE = -2 otherwise). If ITOL is
936 C left unaltered but RTOL and ATOL are uniformly
937 C scaled up by a factor of TOLSF for the next call,
938 C then the solver is deemed likely to succeed.
939 C (The user may also ignore TOLSF and alter the
940 C tolerance parameters in any other way appropriate.)
941 C
942 C NST IWORK(11) the number of steps taken for the problem so far.
943 C
944 C NRE IWORK(12) the number of residual evaluations (RES calls)
945 C for the problem so far.
946 C
947 C NJE IWORK(13) the number of Jacobian evaluations (each involving
948 C an evaluation of a and dr/dy) for the problem so
949 C far. This equals the number of calls to ADDA and
950 C (if MITER = 1) to JAC, and the number of matrix
951 C LU decompositions.
952 C
953 C NQU IWORK(14) the method order last used (successfully).
954 C
955 C NQCUR IWORK(15) the order to be attempted on the next step.
956 C
957 C IMXER IWORK(16) the index of the component of largest magnitude in
958 C the weighted local error vector ( E(i)/EWT(i) ),
959 C on an error return with ISTATE = -4 or -5.
960 C
961 C LENRW IWORK(17) the length of RWORK actually required.
962 C This is defined on normal returns and on an illegal
963 C input return for insufficient storage.
964 C
965 C LENIW IWORK(18) the length of IWORK actually required.
966 C This is defined on normal returns and on an illegal
967 C input return for insufficient storage.
968 C
969 C
970 C The following two arrays are segments of the RWORK array which
971 C may also be of interest to the user as optional outputs.
972 C For each array, the table below gives its internal name,
973 C its base address in RWORK, and its description.
974 C
975 C Name Base Address Description
976 C
977 C YH 21 the Nordsieck history array, of size NYH by
978 C (NQCUR + 1), where NYH is the initial value
979 C of NEQ. For j = 0,1,...,NQCUR, column j+1
980 C of YH contains HCUR**j/factorial(j) times
981 C the j-th derivative of the interpolating
982 C polynomial currently representing the solution,
983 C evaluated at t = TCUR.
984 C
985 C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
986 C corrections on each step, scaled on output to
987 C represent the estimated local error in y on
988 C the last step. This is the vector E in the
989 C description of the error control. It is
990 C defined only on a return from DLSOIBT with
991 C ISTATE = 2.
992 C
993 C-----------------------------------------------------------------------
994 C Part 2. Other Routines Callable.
995 C
996 C The following are optional calls which the user may make to
997 C gain additional capabilities in conjunction with DLSOIBT.
998 C (The routines XSETUN and XSETF are designed to conform to the
999 C SLATEC error handling package.)
1000 C
1001 C Form of Call Function
1002 C CALL XSETUN(LUN) Set the logical unit number, LUN, for
1003 C output of messages from DLSOIBT, if
1004 C the default is not desired.
1005 C The default value of LUN is 6.
1006 C
1007 C CALL XSETF(MFLAG) Set a flag to control the printing of
1008 C messages by DLSOIBT.
1009 C MFLAG = 0 means do not print. (Danger:
1010 C This risks losing valuable information.)
1011 C MFLAG = 1 means print (the default).
1012 C
1013 C Either of the above calls may be made at
1014 C any time and will take effect immediately.
1015 C
1016 C CALL DSRCOM(RSAV,ISAV,JOB) saves and restores the contents of
1017 C the internal Common blocks used by
1018 C DLSOIBT (see Part 3 below).
1019 C RSAV must be a real array of length 218
1020 C or more, and ISAV must be an integer
1021 C array of length 37 or more.
1022 C JOB=1 means save Common into RSAV/ISAV.
1023 C JOB=2 means restore Common from RSAV/ISAV.
1024 C DSRCOM is useful if one is
1025 C interrupting a run and restarting
1026 C later, or alternating between two or
1027 C more problems solved with DLSOIBT.
1028 C
1029 C CALL DINTDY(,,,,,) Provide derivatives of y, of various
1030 C (see below) orders, at a specified point t, if
1031 C desired. It may be called only after
1032 C a successful return from DLSOIBT.
1033 C
1034 C The detailed instructions for using DINTDY are as follows.
1035 C The form of the call is:
1036 C
1037 C CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
1038 C
1039 C The input parameters are:
1040 C
1041 C T = value of independent variable where answers are desired
1042 C (normally the same as the t last returned by DLSOIBT).
1043 C For valid results, T must lie between TCUR - HU and TCUR.
1044 C (See optional outputs for TCUR and HU.)
1045 C K = integer order of the derivative desired. K must satisfy
1046 C 0 .le. K .le. NQCUR, where NQCUR is the current order
1047 C (see optional outputs). The capability corresponding
1048 C to K = 0, i.e. computing y(t), is already provided
1049 C by DLSOIBT directly. Since NQCUR .ge. 1, the first
1050 C derivative dy/dt is always available with DINTDY.
1051 C RWORK(21) = the base address of the history array YH.
1052 C NYH = column length of YH, equal to the initial value of NEQ.
1053 C
1054 C The output parameters are:
1055 C
1056 C DKY = a real array of length NEQ containing the computed value
1057 C of the K-th derivative of y(t).
1058 C IFLAG = integer flag, returned as 0 if K and T were legal,
1059 C -1 if K was illegal, and -2 if T was illegal.
1060 C On an error return, a message is also written.
1061 C-----------------------------------------------------------------------
1062 C Part 3. Common Blocks.
1063 C
1064 C If DLSOIBT is to be used in an overlay situation, the user
1065 C must declare, in the primary overlay, the variables in:
1066 C (1) the call sequence to DLSOIBT, and
1067 C (2) the internal Common block
1068 C /DLS001/ of length 255 (218 double precision words
1069 C followed by 37 integer words),
1070 C
1071 C If DLSOIBT is used on a system in which the contents of internal
1072 C Common blocks are not preserved between calls, the user should
1073 C declare the above Common block in the calling program to insure
1074 C that their contents are preserved.
1075 C
1076 C If the solution of a given problem by DLSOIBT is to be interrupted
1077 C and then later continued, such as when restarting an interrupted run
1078 C or alternating between two or more problems, the user should save,
1079 C following the return from the last DLSOIBT call prior to the
1080 C interruption, the contents of the call sequence variables and the
1081 C internal Common blocks, and later restore these values before the
1082 C next DLSOIBT call for that problem. To save and restore the Common
1083 C blocks, use Subroutine DSRCOM (see Part 2 above).
1084 C
1085 C-----------------------------------------------------------------------
1086 C Part 4. Optionally Replaceable Solver Routines.
1087 C
1088 C Below are descriptions of two routines in the DLSOIBT package which
1089 C relate to the measurement of errors. Either routine can be
1090 C replaced by a user-supplied version, if desired. However, since such
1091 C a replacement may have a major impact on performance, it should be
1092 C done only when absolutely necessary, and only with great caution.
1093 C (Note: The means by which the package version of a routine is
1094 C superseded by the user's version may be system-dependent.)
1095 C
1096 C (a) DEWSET.
1097 C The following subroutine is called just before each internal
1098 C integration step, and sets the array of error weights, EWT, as
1099 C described under ITOL/RTOL/ATOL above:
1100 C SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
1101 C where NEQ, ITOL, RTOL, and ATOL are as in the DLSOIBT call sequence,
1102 C YCUR contains the current dependent variable vector, and
1103 C EWT is the array of weights set by DEWSET.
1104 C
1105 C If the user supplies this subroutine, it must return in EWT(i)
1106 C (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1107 C in y(i) to. The EWT array returned by DEWSET is passed to the DVNORM
1108 C routine (see below), and also used by DLSOIBT in the computation
1109 C of the optional output IMXER, the diagonal Jacobian approximation,
1110 C and the increments for difference quotient Jacobians.
1111 C
1112 C In the user-supplied version of DEWSET, it may be desirable to use
1113 C the current values of derivatives of y. Derivatives up to order NQ
1114 C are available from the history array YH, described above under
1115 C optional outputs. In DEWSET, YH is identical to the YCUR array,
1116 C extended to NQ + 1 columns with a column length of NYH and scale
1117 C factors of H**j/factorial(j). On the first call for the problem,
1118 C given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1119 C NYH is the initial value of NEQ. The quantities NQ, H, and NST
1120 C can be obtained by including in DEWSET the statements:
1121 C DOUBLE PRECISION RLS
1122 C COMMON /DLS001/ RLS(218),ILS(37)
1123 C NQ = ILS(33)
1124 C NST = ILS(34)
1125 C H = RLS(212)
1126 C Thus, for example, the current value of dy/dt can be obtained as
1127 C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is
1128 C unnecessary when NST = 0).
1129 C
1130 C (b) DVNORM.
1131 C The following is a real function routine which computes the weighted
1132 C root-mean-square norm of a vector v:
1133 C D = DVNORM (N, V, W)
1134 C where:
1135 C N = the length of the vector,
1136 C V = real array of length N containing the vector,
1137 C W = real array of length N containing weights,
1138 C D = SQRT( (1/N) * sum(V(i)*W(i))**2 ).
1139 C DVNORM is called with N = NEQ and with W(i) = 1.0/EWT(i), where
1140 C EWT is as set by Subroutine DEWSET.
1141 C
1142 C If the user supplies this function, it should return a non-negative
1143 C value of DVNORM suitable for use in the error control in DLSOIBT.
1144 C None of the arguments should be altered by DVNORM.
1145 C For example, a user-supplied DVNORM routine might:
1146 C -substitute a max-norm of (V(i)*W(i)) for the RMS-norm, or
1147 C -ignore some components of V in the norm, with the effect of
1148 C suppressing the error control on those components of y.
1149 C-----------------------------------------------------------------------
1150 C
1151 C***REVISION HISTORY (YYYYMMDD)
1152 C 19840625 DATE WRITTEN
1153 C 19870330 Major update: corrected comments throughout;
1154 C removed TRET from Common; rewrote EWSET with 4 loops;
1155 C fixed t test in INTDY; added Cray directives in STODI;
1156 C in STODI, fixed DELP init. and logic around PJAC call;
1157 C combined routines to save/restore Common;
1158 C passed LEVEL = 0 in error message calls (except run abort).
1159 C 20010425 Major update: convert source lines to upper case;
1160 C added *DECK lines; changed from 1 to * in dummy dimensions;
1161 C changed names R1MACH/D1MACH to RUMACH/DUMACH;
1162 C renamed routines for uniqueness across single/double prec.;
1163 C converted intrinsic names to generic form;
1164 C removed ILLIN and NTREP (data loaded) from Common;
1165 C removed all 'own' variables from Common;
1166 C changed error messages to quoted strings;
1167 C replaced XERRWV/XERRWD with 1993 revised version;
1168 C converted prologues, comments, error messages to mixed case;
1169 C converted arithmetic IF statements to logical IF statements;
1170 C numerous corrections to prologues and internal comments.
1171 C 20010507 Converted single precision source to double precision.
1172 C 20020502 Corrected declarations in descriptions of user routines.
1173 C 20031105 Restored 'own' variables to Common block, to enable
1174 C interrupt/restart feature.
1175 C 20031112 Added SAVE statements for data-loaded constants.
1176 C 20031117 Changed internal names NRE, LSAVR to NFE, LSAVF resp.
1177 C
1178 C-----------------------------------------------------------------------
1179 C Other routines in the DLSOIBT package.
1180 C
1181 C In addition to Subroutine DLSOIBT, the DLSOIBT package includes the
1182 C following subroutines and function routines:
1183 C DAIGBT computes the initial value of the vector
1184 C dy/dt = A-inverse * g
1185 C DINTDY computes an interpolated value of the y vector at t = TOUT.
1186 C DSTODI is the core integrator, which does one step of the
1187 C integration and the associated error control.
1188 C DCFODE sets all method coefficients and test constants.
1189 C DEWSET sets the error weight vector EWT before each step.
1190 C DVNORM computes the weighted RMS-norm of a vector.
1191 C DSRCOM is a user-callable routine to save and restore
1192 C the contents of the internal Common blocks.
1193 C DPJIBT computes and preprocesses the Jacobian matrix
1194 C and the Newton iteration matrix P.
1195 C DSLSBT manages solution of linear system in chord iteration.
1196 C DDECBT and DSOLBT are routines for solving block-tridiagonal
1197 C systems of linear algebraic equations.
1198 C DGEFA and DGESL are routines from LINPACK for solving full
1199 C systems of linear algebraic equations.
1200 C DDOT is one of the basic linear algebra modules (BLAS).
1201 C DUMACH computes the unit roundoff in a machine-independent manner.
1202 C XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all
1203 C error messages and warnings. XERRWD is machine-dependent.
1204 C Note: DVNORM, DDOT, DUMACH, IXSAV, and IUMACH are function routines.
1205 C All the others are subroutines.
1206 C
1207 C-----------------------------------------------------------------------
1221 LOGICAL IHIT
1222 CHARACTER*60 MSG
1224 C-----------------------------------------------------------------------
1225 C The following internal Common block contains
1226 C (a) variables which are local to any subroutine but whose values must
1227 C be preserved between calls to the routine ("own" variables), and
1228 C (b) variables which are communicated between subroutines.
1229 C The block DLS001 is declared in subroutines DLSOIBT, DINTDY, DSTODI,
1230 C DPJIBT, and DSLSBT.
1231 C Groups of variables are replaced by dummy arrays in the Common
1232 C declarations in routines where those variables are not used.
1233 C-----------------------------------------------------------------------
1240 C
1242 C-----------------------------------------------------------------------
1243 C Block A.
1244 C This code block is executed on every call.
1245 C It tests ISTATE and ITASK for legality and branches appropriately.
1246 C If ISTATE .gt. 1 but the flag INIT shows that initialization has
1247 C not yet been done, an error return occurs.
1248 C If ISTATE = 0 or 1 and TOUT = T, return immediately.
1249 C-----------------------------------------------------------------------
1258 C-----------------------------------------------------------------------
1259 C Block B.
1260 C The next code block is executed for the initial call (ISTATE = 0 or 1)
1261 C or for a continuation call with parameter changes (ISTATE = 3).
1262 C It contains checking of all inputs and various initializations.
1263 C
1264 C First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1265 C MF, MB, and NB.
1266 C-----------------------------------------------------------------------
1282 C Next process and check the optional inputs. --------------------------
1285 MXSTEP = MXSTP0
1286 MXHNIL = MXHNL0
1310 C-----------------------------------------------------------------------
1311 C Set work array pointers and check lengths LRW and LIW.
1312 C Pointers to segments of RWORK and IWORK are named by prefixing L to
1313 C the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1314 C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVR, ACOR.
1315 C-----------------------------------------------------------------------
1330 C Check RTOL and ATOL for legality. ------------------------------------
1340 C If ISTATE = 3, set flag to signal parameter changes to DSTODI. -------
1343 C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into YDOTI.---------
1346 C Reload WM(1) = RWORK(lWM), since lWM may have changed. ---------------
1349 C NEQ was reduced. Zero part of YH to avoid undefined references. -----
1356 C-----------------------------------------------------------------------
1357 C Block C.
1358 C The next block is for the initial call only (ISTATE = 0 or 1).
1359 C It contains all remaining initializations, the call to DAIGBT
1360 C (if ISTATE = 1), and the calculation of the initial step size.
1361 C The error weights in EWT are inverted after being loaded.
1362 C-----------------------------------------------------------------------
1364 TN = T
1383 C Compute initial dy/dt, if necessary, and load it and initial Y into YH
1387 C DLSOIBT must compute initial dy/dt (LYD0 points to YH(*,2)). ---------
1396 C Initial dy/dt was supplied. Load into YH (LYD0 points to YH(*,2).). -
1400 C Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1408 C-----------------------------------------------------------------------
1409 C The coding below computes the step size, H0, to be attempted on the
1410 C first step, unless the user has supplied a value for this.
1411 C First check that TOUT - T differs significantly from zero.
1412 C A scalar tolerance quantity TOL is computed, as MAX(RTOL(i))
1413 C if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted
1414 C so as to be between 100*UROUND and 1.0E-3.
1415 C Then the computed value H0 is given by..
1416 C NEQ
1417 C H0**2 = TOL / ( w0**-2 + (1/NEQ) * Sum ( YDOT(i)/ywt(i) )**2 )
1418 C 1
1419 C where w0 = MAX ( ABS(T), ABS(TOUT) ),
1420 C YDOT(i) = i-th component of initial value of dy/dt,
1421 C ywt(i) = EWT(i)/TOL (a weight for y(i)).
1422 C The sign of H0 is inferred from the initial values of TOUT and T.
1423 C-----------------------------------------------------------------------
1446 C Adjust H0 if necessary to meet HMAX bound. ---------------------------
1449 C Load H with H0 and scale YH(*,2) by H0. ------------------------------
1450 H = H0
1454 C-----------------------------------------------------------------------
1455 C Block D.
1456 C The next code block is for continuation calls only (ISTATE = 2 or 3)
1457 C and is to check stop conditions before taking a step.
1458 C-----------------------------------------------------------------------
1464 T = TOUT
1476 T = TOUT
1487 C-----------------------------------------------------------------------
1488 C Block E.
1489 C The next block is normally executed for all calls and contains
1490 C the call to the one-step core integrator DSTODI.
1491 C
1492 C This is a looping point for the integration steps.
1493 C
1494 C First check for too many steps being taken, update EWT (if not at
1495 C start of problem), check for too much accuracy being requested, and
1496 C check for H below the roundoff level in T.
1497 C-----------------------------------------------------------------------
1524 C-----------------------------------------------------------------------
1525 C CALL DSTODI(NEQ,Y,YH,NYH,YH1,EWT,SAVF,SAVR,ACOR,WM,IWM,RES,
1526 C ADDA,JAC,DPJIBT,DSLSBT)
1527 C Note: SAVF in DSTODI occupies the same space as YDOTI in DLSOIBT.
1528 C-----------------------------------------------------------------------
1534 C
1535 C KGO = 1:success; 2:error test failure; 3:convergence failure;
1536 C 4:RES ordered return; 5:RES returned error.
1537 C-----------------------------------------------------------------------
1538 C Block F.
1539 C The following block handles the case of a successful return from the
1540 C core integrator (KFLAG = 0). Test for stop conditions.
1541 C-----------------------------------------------------------------------
1544 C ITASK = 1. If TOUT has been reached, interpolate. -------------------
1547 T = TOUT
1549 C ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1552 C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1555 T = TOUT
1565 C ITASK = 5. see if TCRIT was reached and jump to exit. ---------------
1568 C-----------------------------------------------------------------------
1569 C Block G.
1570 C The following block handles all successful returns from DLSOIBT.
1571 C If ITASK .ne. 1, Y is loaded from YH and T is set accordingly.
1572 C ISTATE is set to 2, and the optional outputs are loaded into the
1573 C work arrays before returning.
1574 C-----------------------------------------------------------------------
1577 T = TN
1590 RETURN
1591 C-----------------------------------------------------------------------
1592 C Block H.
1593 C The following block handles all unsuccessful returns other than
1594 C those for illegal input. First the error message routine is called.
1595 C If there was an error test or convergence test failure, IMXER is set.
1596 C Then Y is loaded from YH and T is set to TN.
1597 C The optional outputs are loaded into the work arrays before returning.
1598 C-----------------------------------------------------------------------
1599 C The maximum number of steps was taken before reaching TOUT. ----------
1606 C EWT(i) .le. 0.0 for some i (not at start of problem). ----------------
1612 C Too much accuracy requested for machine precision. -------------------
1620 C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1627 C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1636 C IRES = 3 returned by RES, despite retries by DSTODI.------------------
1643 C DAIGBT failed because a diagonal block of A matrix was singular. -----
1650 RETURN
1651 C DAIGBT failed because RES set IRES to 2 or 3. ------------------------
1659 RETURN
1660 C Compute IMXER if relevant. -------------------------------------------
1666 BIG = SIZE
1667 IMXER = I
1670 C Compute residual if relevant. ----------------------------------------
1684 C Set Y vector, T, and optional outputs. -------------------------------
1696 RETURN
1697 C-----------------------------------------------------------------------
1698 C Block I.
1699 C The following block handles all error returns due to illegal input
1700 C (ISTATE = -3), as detected before calling the core integrator.
1701 C First the error message routine is called. If the illegal input
1702 C is a negative ISTATE, the run is aborted (apparent infinite loop).
1703 C-----------------------------------------------------------------------
1791 C
1793 RETURN
1794 C
1797 RETURN
1798 C----------------------- End of Subroutine DLSOIBT ---------------------
1799 END