1
subroutine lsoibt (res, adda, jac, neq, y, ydoti, t, tout, itol,
2
1 rtol, atol, itask, istate, iopt, rwork, lrw, iwork, liw, mf )
3
external res, adda, jac
4
integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
5
double precision y, ydoti, t, tout, rtol, atol, rwork
6
dimension neq(1), y(1), ydoti(1), rtol(1), atol(1), rwork(lrw),
8
c-----------------------------------------------------------------------
9
c this is the march 30, 1987 version of lsoibt..
10
c livermore solver for ordinary differential equations given in
11
c implicit form, with block-tridiagonal jacobian treatment.
12
c this version is in double precision.
14
c lsoibt solves the initial value problem for linearly implicit
15
c systems of first order ode-s,
16
c a(t,y) * dy/dt = g(t,y) , where a(t,y) is a square matrix,
17
c or, in component form,
18
c ( a * ( dy / dt )) + ... + ( a * ( dy / dt )) =
21
c = g ( t, y , y ,..., y ) ( i = 1,...,neq )
24
c if a is singular, this is a differential-algebraic system.
26
c lsoibt is a variant version of the lsodi package, for the case where
27
c the matrices a, dg/dy, and d(a*s)/dy are all block-tridiagonal.
28
c-----------------------------------------------------------------------
30
c alan c. hindmarsh, odepack, a systematized collection of ode
31
c solvers, in scientific computing, r. s. stepleman et al. (eds.),
32
c north-holland, amsterdam, 1983, pp. 55-64.
33
c-----------------------------------------------------------------------
34
c authors.. charles s. kenney
36
c naval weapons center
37
c china lake, ca 93555
39
c jeffrey f. painter and
41
c computing and mathematics research division, l-316
42
c lawrence livermore national laboratory
43
c livermore, ca 94550.
44
c-----------------------------------------------------------------------
47
c communication between the user and the lsoibt 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.
54
c a. first, provide a subroutine of the form..
55
c subroutine res (neq, t, y, s, r, ires)
56
c dimension y(neq), s(neq), r(neq)
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 lsoibt. (for uses of ires, see the
64
c paragraph on res in the full description below.)
66
c b. next, identify the block structure of the matrices a = a(t,y) and
67
c dr/dy. lsoibt 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 lsoibt 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.)
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 dimension y(neq), 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..
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)
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 dimension y(neq), s(neq), pa(mb,mb,nb), pb(mb,mb,nb), pc(mb,mb,nb)
126
c which computes dr/dy as a function of t, y, and s. here t, y, and
127
c s are inputs, and the routine is to load dr/dy into pa, pb, pc,
128
c according to the structure description in paragraph b above.
129
c that is, load the diagonal blocks into pa, the superdiagonal blocks
130
c (and block (nb,nb-2) ) into pb, and the subdiagonal blocks (and
131
c block (1,3) ) into pc. the blocks in block-row k of dr/dy are to
132
c be loaded into pa(*,*,k), pb(*,*,k), and pc(*,*,k).
133
c only nonzero elements need be loaded, and the indexing
134
c of pa, pb, and pc is the same as in the adda routine.
135
c note that if a is independent of y (or this dependence
136
c is weak enough to be ignored) then jac is to compute dg/dy.
137
c if it is not feasible to provide a jac routine, use
138
c mf = 22, and lsoibt will compute an approximate jacobian
139
c internally by difference quotients.
141
c e. next decide whether or not to provide the initial value of the
142
c derivative vector dy/dt. if the initial value of a(t,y) is
143
c nonsingular (and not too ill-conditioned), you may let lsoibt compute
144
c this vector (istate = 0). (lsoibt will solve the system a*s = g for
145
c s, with initial values of a and g.) if a(t,y) is initially
146
c singular, then the system is a differential-algebraic system, and
147
c you must make use of the particular form of the system to compute the
148
c initial values of y and dy/dt. in that case, use istate = 1 and
149
c load the initial value of dy/dt into the array ydoti.
150
c the input array ydoti and the initial y array must be consistent with
151
c the equations a*dy/dt = g. this implies that the initial residual
152
c r = g(t,y) - a(t,y)*ydoti must be approximately zero.
154
c f. write a main program which calls subroutine lsoibt once for
155
c each point at which answers are desired. this should also provide
156
c for possible use of logical unit 6 for output of error messages
157
c by lsoibt. on the first call to lsoibt, supply arguments as follows..
158
c res = name of user subroutine for residual function r.
159
c adda = name of user subroutine for computing and adding a(t,y).
160
c jac = name of user subroutine for jacobian matrix dr/dy
161
c (mf = 21). if not used, pass a dummy name.
162
c note.. the names for the res and adda routines and (if used) the
163
c jac routine must be declared external in the calling program.
164
c neq = number of scalar equations in the system.
165
c y = array of initial values, of length neq.
166
c ydoti = array of length neq (containing initial dy/dt if istate = 1).
167
c t = the initial value of the independent variable.
168
c tout = first point where output is desired (.ne. t).
169
c itol = 1 or 2 according as atol (below) is a scalar or array.
170
c rtol = relative tolerance parameter (scalar).
171
c atol = absolute tolerance parameter (scalar or array).
172
c the estimated local error in y(i) will be controlled so as
173
c to be roughly less (in magnitude) than
174
c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or
175
c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2.
176
c thus the local error test passes if, in each component,
177
c either the absolute error is less than atol (or atol(i)),
178
c or the relative error is less than rtol.
179
c use rtol = 0.0 for pure absolute error control, and
180
c use atol = 0.0 (or atol(i) = 0.0) for pure relative error
181
c control. caution.. actual (global) errors may exceed these
182
c local tolerances, so choose them conservatively.
183
c itask = 1 for normal computation of output values of y at t = tout.
184
c istate = integer flag (input and output). set istate = 1 if the
185
c initial dy/dt is supplied, and 0 otherwise.
186
c iopt = 0 to indicate no optional inputs used.
187
c rwork = real work array of length at least..
188
c 22 + 9*neq + 3*mb*mb*nb for mf = 21 or 22.
189
c lrw = declared length of rwork (in user-s dimension).
190
c iwork = integer work array of length at least 20 + neq.
191
c input in iwork(1) the block size mb and in iwork(2) the
192
c number nb of blocks in each direction along the matrix a.
193
c these must satisfy mb .ge. 1, nb .ge. 4, and mb*nb = neq.
194
c liw = declared length of iwork (in user-s dimension).
195
c mf = method flag. standard values are..
196
c 21 for a user-supplied jacobian.
197
c 22 for an internally generated jacobian.
198
c for other choices of mf, see the paragraph on mf in
199
c the full description below.
200
c note that the main program must declare arrays y, ydoti, rwork, iwork,
203
c g. the output from the first call (or any call) is..
204
c y = array of computed values of y(t) vector.
205
c t = corresponding value of independent variable (normally tout).
206
c istate = 2 if lsoibt was successful, negative otherwise.
207
c -1 means excess work done on this call (check all inputs).
208
c -2 means excess accuracy requested (tolerances too small).
209
c -3 means illegal input detected (see printed message).
210
c -4 means repeated error test failures (check all inputs).
211
c -5 means repeated convergence failures (perhaps bad jacobian
212
c supplied or wrong choice of tolerances).
213
c -6 means error weight became zero during problem. (solution
214
c component i vanished, and atol or atol(i) = 0.)
215
c -7 cannot occur in casual use.
216
c -8 means lsoibt was unable to compute the initial dy/dt.
217
c in casual use, this means a(t,y) is initially singular.
218
c supply ydoti and use istate = 1 on the first call.
220
c if lsoibt returns istate = -1, -4, or -5, then the output of
221
c lsoibt also includes ydoti = array containing residual vector
222
c r = g - a * dy/dt evaluated at the current t, y, and dy/dt.
224
c h. to continue the integration after a successful return, simply
225
c reset tout and call lsoibt again. no other parameters need be reset.
228
c-----------------------------------------------------------------------
231
c the following is an example problem, with the coding needed
232
c for its solution by lsoibt. 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 p.d.e. 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.
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)
280
c call lsoibt (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(7h at t =,f5.2,14h no. steps =,i4,11h no. r-s =,i4,
284
c 1 11h no. j-s =,i3)
285
c if (istate .ne. 2) go to 90
286
c tout = tout + 0.1d0
288
c write(6,50) (y(i),i=1,neq)
289
c 50 format(/24h final solution values../9(5e12.4/))
291
c 90 write(6,95) istate
292
c 95 format(///22h error halt.. istate =,i3)
296
c subroutine resid (n, t, y, s, r, ires)
297
c double precision delx, eta, eodsq, r, s, t, y
298
c dimension y(n), s(n), r(n)
299
c data eta/0.05d0/, delx/0.05d0/
300
c eodsq = eta/delx**2
301
c r(1) = eodsq*(y(3) - 2.0d0*y(2) + y(1)) - s(1)
304
c r(i) = (y(i-1)**2 - y(i+1)**2)/(4.0d0*delx)
305
c 1 + eodsq*(y(i+1) - 2.0d0*y(i) + y(i-1))
306
c 2 - (s(i-1) + 4.0d0*s(i) + s(i+1))/6.0d0
308
c r(n) = eodsq*(y(n-2) - 2.0d0*y(nm1) + y(n)) - s(n)
312
c subroutine addabt (n, t, y, mb, nb, pa, pb, pc)
313
c double precision pa, pb, pc, t, y
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.0d0
318
c pa(1,1,k) = pa(1,1,k) + (4.0d0/6.0d0)
319
c pb(1,1,k) = pb(1,1,k) + (1.0d0/6.0d0)
320
c pc(1,1,k) = pc(1,1,k) + (1.0d0/6.0d0)
322
c pa(1,1,n) = pa(1,1,n) + 1.0d0
326
c subroutine jacbt (n, t, y, s, mb, nb, pa, pb, pc)
327
c double precision delx, eta, eodsq, pa, pb, pc, s, t, y
328
c dimension y(n), s(n), pa(mb,mb,nb),pb(mb,mb,nb),pc(mb,mb,nb)
329
c data eta/0.05d0/, delx/0.05d0/
330
c eodsq = eta/delx**2
332
c pb(1,1,1) = -2.0d0*eodsq
335
c pa(1,1,k) = -2.0d0*eodsq
336
c pb(1,1,k) = -y(k+1)*(0.5d0/delx) + eodsq
337
c pc(1,1,k) = y(k-1)*(0.5d0/delx) + eodsq
340
c pc(1,1,n) = -2.0d0*eodsq
345
c the output of this program (on a cdc-7600 in single precision)
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
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
363
c-----------------------------------------------------------------------
364
c full description of user interface to lsoibt.
366
c the user interface to lsoibt consists of the following parts.
368
c i. the call sequence to subroutine lsoibt, which is a driver
369
c routine for the solver. this includes descriptions of both
370
c the call sequence arguments and of user-supplied routines.
371
c following these descriptions is a description of
372
c optional inputs available through the call sequence, and then
373
c a description of optional outputs (in the work arrays).
375
c ii. descriptions of other routines in the lsoibt package that may be
376
c (optionally) called by the user. these provide the ability to
377
c alter error message handling, save and restore the internal
378
c common, and obtain specified derivatives of the solution y(t).
380
c iii. descriptions of common blocks to be declared in overlay
381
c or similar environments, or to be saved when doing an interrupt
382
c of the problem and continued solution later.
384
c iv. description of two routines in the lsoibt package, either of
385
c which the user may replace with his own version, if desired.
386
c these relate to the measurement of errors.
388
c-----------------------------------------------------------------------
389
c part i. call sequence.
391
c the call sequence parameters used for input only are
392
c res, adda, jac, neq, tout, itol, rtol, atol, itask,
393
c iopt, lrw, liw, mf,
394
c and those used for both input and output are
395
c y, t, istate, ydoti.
396
c the work arrays rwork and iwork are also used for additional and
397
c optional inputs and optional outputs. (the term output here refers
398
c to the return from subroutine lsoibt to the user-s calling program.)
400
c the legality of input parameters will be thoroughly checked on the
401
c initial call for the problem, but not checked thereafter unless a
402
c change in input parameters is flagged by istate = 3 on input.
404
c the descriptions of the call arguments are as follows.
406
c res = the name of the user-supplied subroutine which supplies
407
c the residual vector for the ode system, defined by
408
c r = g(t,y) - a(t,y) * s
409
c as a function of the scalar t and the vectors
410
c s and y ( s approximates dy/dt ). this
411
c subroutine is to have the form
412
c subroutine res ( neq, t, y, s, r, ires )
413
c dimension y(1), s(1), r(1)
414
c where neq, t, y, s, and ires are input, and r and
415
c ires are output. y, s, and r are arrays of length neq.
416
c in dimension statements such as that above, 1 is a
417
c dummy dimension. it can be replaced by any value.
418
c on input, ires indicates how lsoibt will use the
419
c returned array r, as follows..
420
c ires = 1 means that lsoibt needs the full residual,
421
c r = g - a*s, exactly.
422
c ires = -1 means that lsoibt is using r only to compute
423
c the jacobian dr/dy by difference quotients.
424
c the res routine can ignore ires, or it can omit some terms
425
c if ires = -1. if a does not depend on y, then res can
426
c just return r = g when ires = -1. if g - a*s contains other
427
c additive terms that are independent of y, these can also be
428
c dropped, if done consistently, when ires = -1.
429
c the subroutine should set the flag ires if it
430
c encounters a halt condition or illegal input.
431
c otherwise, it should not reset ires. on output,
432
c ires = 1 or -1 represents a normal return, and
433
c lsoibt continues integrating the ode. leave ires
434
c unchanged from its input value.
435
c ires = 2 tells lsoibt to immediately return control
436
c to the calling program, with istate = 3. this lets
437
c the calling program change parameters of the prob-
439
c ires = 3 represents an error condition (for example, an
440
c illegal value of y). lsoibt tries to integrate the ode
441
c without getting ires = 3 from res. if it cannot, lsoibt
442
c returns with istate = -7 or -1.
443
c on an lsoibt return with istate = 3, -1, or -7, the values
444
c of t and y returned correspond to the last point reached
445
c successfully without getting the flag ires = 2 or 3.
446
c the flag values ires = 2 and 3 should not be used to
447
c handle switches or root-stop conditions. this is better
448
c done by calling lsoibt in a one-step mode and checking the
449
c stopping function for a sign change at each step.
450
c if quantities computed in the res routine are needed
451
c externally to lsoibt, an extra call to res should be made
452
c for this purpose, for consistent and accurate results.
453
c to get the current dy/dt for the s argument, use intdy.
454
c res must be declared external in the calling
455
c program. see note below for more about res.
457
c adda = the name of the user-supplied subroutine which adds the
458
c matrix a = a(t,y) to another matrix, p, stored in
459
c block-tridiagonal form. this routine is to have the form
460
c subroutine adda (neq, t, y, mb, nb, pa, pb, pc)
462
c 1 pa(mb,mb,nb), pb(mb,mb,nb), pc(mb,mb,nb)
463
c where neq, t, y, mb, nb, and the arrays pa, pb, and pc
464
c are input, and the arrays pa, pb, and pc are output.
465
c y is an array of length neq, and the arrays pa, pb, pc
466
c are all mb by mb by nb.
467
c here a block-tridiagonal structure is assumed for a(t,y),
468
c and also for the matrix p to which a is added here,
469
c as described in paragraph b of the summary of usage above.
470
c thus the affect of adda should be the following..
474
c pa(i,j,k) = pa(i,j,k) +
475
c ( (i,j) element of k-th diagonal block of a)
476
c pb(i,j,k) = pb(i,j,k) +
477
c ( (i,j) element of block (k,k+1) of a,
478
c or block (nb,nb-2) if k = nb)
479
c pc(i,j,k) = pc(i,j,k) +
480
c ( (i,j) element of block (k,k-1) of a,
481
c or block (1,3) if k = 1)
485
c adda must be declared external in the calling program.
486
c see note below for more information about adda.
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 -- otherwise a dummy name can be
491
c passed. this subroutine is to have the form
492
c subroutine jac (neq, t, y, s, mb, nb, pa, pb, pc)
493
c dimension y(neq), s(neq),
494
c 1 pa(mb,mb,nb), pb(mb,mb,nb), pc(mb,mb,nb)
495
c where neq, t, y, s, mb, nb, and the arrays pa, pb, and pc
496
c are input, and the arrays pa, pb, and pc are output.
497
c y and s are arrays of length neq, and the arrays pa, pb, pc
498
c are all mb by mb by nb.
499
c pa, pb, and pc are to be loaded with partial derivatives
500
c (elements of the jacobian matrix) on output, in terms of the
501
c block-tridiagonal structure assumed, as described
502
c in paragraph b of the summary of usage above.
503
c that is, load the diagonal blocks into pa, the
504
c superdiagonal blocks (and block (nb,nb-2) ) intp pb, and
505
c the subdiagonal blocks (and block (1,3) ) into pc.
506
c the blocks in block-row k of dr/dy are to be loaded into
507
c pa(*,*,k), pb(*,*,k), and pc(*,*,k).
508
c thus the affect of jac should be the following..
512
c pa(i,j,k) = ( (i,j) element of
513
c k-th diagonal block of dr/dy)
514
c pb(i,j,k) = ( (i,j) element of block (k,k+1)
515
c of dr/dy, or block (nb,nb-2) if k = nb)
516
c pc(i,j,k) = ( (i,j) element of block (k,k-1)
517
c of dr/dy, or block (1,3) if k = 1)
521
c pa, pb, and pc are preset to zero by the solver,
522
c so that only the nonzero elements need be loaded by jac.
523
c each call to jac is preceded by a call to res with the same
524
c arguments neq, t, y, and s. thus to gain some efficiency,
525
c intermediate quantities shared by both calculations may be
526
c saved in a user common block by res and not recomputed by jac
527
c if desired. also, jac may alter the y array, if desired.
528
c jac need not provide dr/dy exactly. a crude
529
c approximation will do, so that lsoibt may be used when
530
c a and dr/dy are not really block-tridiagonal, but are close
531
c to matrices that are.
532
c jac must be declared external in the calling program.
533
c see note below for more about jac.
535
c note on res, adda, and jac-- these
536
c subroutines may access user-defined quantities in
537
c neq(2),... and/or in y(neq(1)+1),... if neq is an array
538
c (dimensioned in the subroutines) and/or y has length
539
c exceeding neq(1). however, these routines should not alter
540
c neq(1), y(1),...,y(neq) or any other input variables.
541
c see the descriptions of neq and y below.
543
c neq = the size of the system (number of first order ordinary
544
c differential equations or scalar algebraic equations).
545
c used only for input.
546
c neq may be decreased, but not increased, during the problem.
547
c if neq is decreased (with istate = 3 on input), the
548
c remaining components of y should be left undisturbed, if
549
c these are to be accessed in res, adda, or jac.
551
c normally, neq is a scalar, and it is generally referred to
552
c as a scalar in this user interface description. however,
553
c neq may be an array, with neq(1) set to the system size.
554
c (the lsoibt package accesses only neq(1).) in either case,
555
c this parameter is passed as the neq argument in all calls
556
c to res, adda, and jac. hence, if it is an array,
557
c locations neq(2),... may be used to store other integer data
558
c and pass it to res, adda, or jac. each such subroutine
559
c must include neq in a dimension statement in that case.
561
c y = a real array for the vector of dependent variables, of
562
c length neq or more. used for both input and output on the
563
c first call (istate = 0 or 1), and only for output on other
564
c calls. on the first call, y must contain the vector of
565
c initial values. on output, y contains the computed solution
566
c vector, evaluated at t. if desired, the y array may be used
567
c for other purposes between calls to the solver.
569
c this array is passed as the y argument in all calls to res,
570
c adda, and jac. hence its length may exceed neq,
571
c and locations y(neq+1),... may be used to store other real
572
c data and pass it to res, adda, or jac. (the lsoibt
573
c package accesses only y(1),...,y(neq). )
575
c ydoti = a real array for the initial value of the vector
576
c dy/dt and for work space, of dimension at least neq.
579
c if istate = 0 then lsoibt will compute the initial value
580
c of dy/dt, if a is nonsingular. thus ydoti will
581
c serve only as work space and may have any value.
582
c if istate = 1 then ydoti must contain the initial value
584
c if istate = 2 or 3 (continuation calls) then ydoti
585
c may have any value.
586
c n.b.- if the initial value of a is singular, then
587
c lsoibt cannot compute the initial value of dy/dt, so
588
c it must be provided in ydoti, with istate=1.
590
c on output, when lsoibt terminates abnormally with istate =
591
c -1, -4, or -5, ydoti will contain the residual
592
c r = g(t,y) - a(t,y)*(dy/dt). if r is large, t is near
593
c its initial value, and ydoti is supplied with istate=1,
594
c there may have been an incorrect input value of
595
c ydoti = dy/dt or the problem ( as given to lsoibt )
596
c may not have a solution.
598
c if desired, the ydoti array may be used for other
599
c purposes between calls to the solver.
601
c t = the independent variable. on input, t is used only on the
602
c first call, as the initial point of the integration.
603
c on output, after each call, t is the value at which a
604
c computed solution y is evaluated (usually the same as tout).
605
c on an error return, t is the farthest point reached.
607
c tout = the next value of t at which a computed solution is desired.
608
c used only for input.
610
c when starting the problem (istate = 0 or 1), tout may be
611
c equal to t for one call, then should .ne. t for the next
612
c call. for the initial t, an input value of tout .ne. t is
613
c used in order to determine the direction of the integration
614
c (i.e. the algebraic sign of the step sizes) and the rough
615
c scale of the problem. integration in either direction
616
c (forward or backward in t) is permitted.
618
c if itask = 2 or 5 (one-step modes), tout is ignored after
619
c the first call (i.e. the first call with tout .ne. t).
620
c otherwise, tout is required on every call.
622
c if itask = 1, 3, or 4, the values of tout need not be
623
c monotone, but a value of tout which backs up is limited
624
c to the current internal t interval, whose endpoints are
625
c tcur - hu and tcur (see optional outputs, below, for
628
c itol = an indicator for the type of error control. see
629
c description below under atol. used only for input.
631
c rtol = a relative error tolerance parameter, either a scalar or
632
c an array of length neq. see description below under atol.
635
c atol = an absolute error tolerance parameter, either a scalar or
636
c an array of length neq. input only.
638
c the input parameters itol, rtol, and atol determine
639
c the error control performed by the solver. the solver will
640
c control the vector e = (e(i)) of estimated local errors
641
c in y, according to an inequality of the form
642
c rms-norm of ( e(i)/ewt(i) ) .le. 1,
643
c where ewt(i) = rtol(i)*abs(y(i)) + atol(i),
644
c and the rms-norm (root-mean-square norm) here is
645
c rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i))
646
c is a vector of weights which must always be positive, and
647
c the values of rtol and atol should all be non-negative.
648
c the following table gives the types (scalar/array) of
649
c rtol and atol, and the corresponding form of ewt(i).
651
c itol rtol atol ewt(i)
652
c 1 scalar scalar rtol*abs(y(i)) + atol
653
c 2 scalar array rtol*abs(y(i)) + atol(i)
654
c 3 array scalar rtol(i)*abs(y(i)) + atol
655
c 4 array scalar rtol(i)*abs(y(i)) + atol(i)
657
c when either of these parameters is a scalar, it need not
658
c be dimensioned in the user-s calling program.
660
c if none of the above choices (with itol, rtol, and atol
661
c fixed throughout the problem) is suitable, more general
662
c error controls can be obtained by substituting
663
c user-supplied routines for the setting of ewt and/or for
664
c the norm calculation. see part iv below.
666
c if global errors are to be estimated by making a repeated
667
c run on the same problem with smaller tolerances, then all
668
c components of rtol and atol (i.e. of ewt) should be scaled
671
c itask = an index specifying the task to be performed.
672
c input only. itask has the following values and meanings.
673
c 1 means normal computation of output values of y(t) at
674
c t = tout (by overshooting and interpolating).
675
c 2 means take one step only and return.
676
c 3 means stop at the first internal mesh point at or
677
c beyond t = tout and return.
678
c 4 means normal computation of output values of y(t) at
679
c t = tout but without overshooting t = tcrit.
680
c tcrit must be input as rwork(1). tcrit may be equal to
681
c or beyond tout, but not behind it in the direction of
682
c integration. this option is useful if the problem
683
c has a singularity at or beyond t = tcrit.
684
c 5 means take one step, without passing tcrit, and return.
685
c tcrit must be input as rwork(1).
687
c note.. if itask = 4 or 5 and the solver reaches tcrit
688
c (within roundoff), it will return t = tcrit (exactly) to
689
c indicate this (unless itask = 4 and tout comes before tcrit,
690
c in which case answers at t = tout are returned first).
692
c istate = an index used for input and output to specify the
693
c state of the calculation.
695
c on input, the values of istate are as follows.
696
c 0 means this is the first call for the problem, and
697
c lsoibt is to compute the initial value of dy/dt
698
c (while doing other initializations). see note below.
699
c 1 means this is the first call for the problem, and
700
c the initial value of dy/dt has been supplied in
701
c ydoti (lsoibt will do other initializations). see note
703
c 2 means this is not the first call, and the calculation
704
c is to continue normally, with no change in any input
705
c parameters except possibly tout and itask.
706
c (if itol, rtol, and/or atol are changed between calls
707
c with istate = 2, the new values will be used but not
708
c tested for legality.)
709
c 3 means this is not the first call, and the
710
c calculation is to continue normally, but with
711
c a change in input parameters other than
712
c tout and itask. changes are allowed in
713
c neq, itol, rtol, atol, iopt, lrw, liw, mf, mb, nb,
714
c and any of the optional inputs except h0.
715
c (see iwork description for mb and nb.)
716
c note.. a preliminary call with tout = t is not counted
717
c as a first call here, as no initialization or checking of
718
c input is done. (such a call is sometimes useful for the
719
c purpose of outputting the initial conditions.)
720
c thus the first call for which tout .ne. t requires
721
c istate = 0 or 1 on input.
723
c on output, istate has the following values and meanings.
724
c 0 or 1 means nothing was done, as tout was equal to t with
725
c istate = 0 or 1 on input. (however, an internal counter
726
c was set to detect and prevent repeated calls of this
728
c 2 means that the integration was performed successfully.
729
c 3 means that the user-supplied subroutine res signalled
730
c lsoibt to halt the integration and return (ires=2).
731
c integration as far as t was achieved with no occurrence
732
c of ires=2, but this flag was set on attempting the next
734
c -1 means an excessive amount of work (more than mxstep
735
c steps) was done on this call, before completing the
736
c requested task, but the integration was otherwise
737
c successful as far as t. (mxstep is an optional input
738
c and is normally 500.) to continue, the user may
739
c simply reset istate to a value .gt. 1 and call again
740
c (the excess work step counter will be reset to 0).
741
c in addition, the user may increase mxstep to avoid
742
c this error return (see below on optional inputs).
743
c -2 means too much accuracy was requested for the precision
744
c of the machine being used. this was detected before
745
c completing the requested task, but the integration
746
c was successful as far as t. to continue, the tolerance
747
c parameters must be reset, and istate must be set
748
c to 3. the optional output tolsf may be used for this
749
c purpose. (note.. if this condition is detected before
750
c taking any steps, then an illegal input return
751
c (istate = -3) occurs instead.)
752
c -3 means illegal input was detected, before taking any
753
c integration steps. see written message for details.
754
c note.. if the solver detects an infinite loop of calls
755
c to the solver with illegal input, it will cause
757
c -4 means there were repeated error test failures on
758
c one attempted step, before completing the requested
759
c task, but the integration was successful as far as t.
760
c the problem may have a singularity, or the input
761
c may be inappropriate.
762
c -5 means there were repeated convergence test failures on
763
c one attempted step, before completing the requested
764
c task, but the integration was successful as far as t.
765
c this may be caused by an inaccurate jacobian matrix.
766
c -6 means ewt(i) became zero for some i during the
767
c integration. pure relative error control (atol(i)=0.0)
768
c was requested on a variable which has now vanished.
769
c the integration was successful as far as t.
770
c -7 means that the user-supplied subroutine res set
771
c its error flag (ires=3) despite repeated tries by lsoibt
772
c to avoid that condition.
773
c -8 means that istate was 0 on input but lsoibt was unable
774
c to compute the initial value of dy/dt. see the
775
c printed message for details.
777
c note.. since the normal output value of istate is 2,
778
c it does not need to be reset for normal continuation.
779
c similarly, istate need not be reset if res told lsoibt
780
c to return because the calling program must change
781
c the parameters of the problem.
782
c also, since a negative input value of istate will be
783
c regarded as illegal, a negative output value requires the
784
c user to change it, and possibly other inputs, before
785
c calling the solver again.
787
c iopt = an integer flag to specify whether or not any optional
788
c inputs are being used on this call. input only.
789
c the optional inputs are listed separately below.
790
c iopt = 0 means no optional inputs are being used.
791
c default values will be used in all cases.
792
c iopt = 1 means one or more optional inputs are being used.
794
c rwork = a real working array (double precision).
795
c the length of rwork must be at least
796
c 20 + nyh*(maxord + 1) + 3*neq + lenwm where
797
c nyh = the initial value of neq,
798
c maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
799
c smaller value is given as an optional input),
800
c lenwm = 3*mb*mb*nb + 2.
801
c (see mf description for the definition of meth.)
802
c thus if maxord has its default value and neq is constant,
804
c 22 + 16*neq + 3*mb*mb*nb for mf = 11 or 12,
805
c 22 + 9*neq + 3*mb*mb*nb for mf = 21 or 22.
806
c the first 20 words of rwork are reserved for conditional
807
c and optional inputs and optional outputs.
809
c the following word in rwork is a conditional input..
810
c rwork(1) = tcrit = critical value of t which the solver
811
c is not to overshoot. required if itask is
812
c 4 or 5, and ignored otherwise. (see itask.)
814
c lrw = the length of the array rwork, as declared by the user.
815
c (this will be checked by the solver.)
817
c iwork = an integer work array. the length of iwork must be at least
818
c 20 + neq . the first few words of iwork are used for
819
c additional and optional inputs and optional outputs.
821
c the following 2 words in iwork are additional required
823
c iwork(1) = mb = block size
824
c iwork(2) = nb = number of blocks in the main diagonal
825
c these must satisfy mb .ge. 1, nb .ge. 4, and mb*nb = neq.
827
c liw = the length of the array iwork, as declared by the user.
828
c (this will be checked by the solver.)
830
c note.. the work arrays must not be altered between calls to lsoibt
831
c for the same problem, except possibly for the additional and
832
c optional inputs, and except for the last 3*neq words of rwork.
833
c the latter space is used for internal scratch space, and so is
834
c available for use by the user outside lsoibt between calls, if
835
c desired (but not for use by res, adda, or jac).
837
c mf = the method flag. used only for input. the legal values of
838
c mf are 11, 12, 21, and 22.
839
c mf has decimal digits meth and miter.. mf = 10*meth + miter.
840
c meth indicates the basic linear multistep method..
841
c meth = 1 means the implicit adams method.
842
c meth = 2 means the method based on backward
843
c differentiation formulas (bdf-s).
844
c the bdf method is strongly preferred for stiff prob-
845
c lems, while the adams method is preferred when the prob-
846
c lem is not stiff. if the matrix a(t,y) is nonsingular,
847
c stiffness here can be taken to mean that of the explicit
848
c ode system dy/dt = a**(-1) * g. if a is singular, the
849
c concept of stiffness is not well defined.
850
c if you do not know whether the problem is stiff, we
851
c recommend using meth = 2. if it is stiff, the advan-
852
c tage of meth = 2 over 1 will be great, while if it is
853
c not stiff, the advantage of meth = 1 will be slight.
854
c if maximum efficiency is important, some experimentation
855
c with meth may be necessary.
856
c miter indicates the corrector iteration method..
857
c miter = 1 means chord iteration with a user-supplied
858
c block-tridiagonal jacobian.
859
c miter = 2 means chord iteration with an internally
860
c generated (difference quotient) block-
861
c tridiagonal jacobian approximation, using
862
c 3*mb+1 extra calls to res per dr/dy evaluation.
863
c if miter = 1, the user must supply a subroutine
864
c jac (the name is arbitrary) as described above under jac.
865
c for miter = 2, a dummy argument can be used.
866
c-----------------------------------------------------------------------
869
c the following is a list of the optional inputs provided for in the
870
c call sequence. (see also part ii.) for each such input variable,
871
c this table lists its name as used in this documentation, its
872
c location in the call sequence, its meaning, and the default value.
873
c the use of any of these inputs requires iopt = 1, and in that
874
c case all of these inputs are examined. a value of zero for any
875
c of these optional inputs will cause the default value to be used.
876
c thus to use a subset of the optional inputs, simply preload
877
c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and
878
c then set those of interest to nonzero values.
880
c name location meaning and default value
882
c h0 rwork(5) the step size to be attempted on the first step.
883
c the default value is determined by the solver.
885
c hmax rwork(6) the maximum absolute step size allowed.
886
c the default value is infinite.
888
c hmin rwork(7) the minimum absolute step size allowed.
889
c the default value is 0. (this lower bound is not
890
c enforced on the final step before reaching tcrit
891
c when itask = 4 or 5.)
893
c maxord iwork(5) the maximum order to be allowed. the default
894
c value is 12 if meth = 1, and 5 if meth = 2.
895
c if maxord exceeds the default value, it will
896
c be reduced to the default value.
897
c if maxord is changed during the problem, it may
898
c cause the current order to be reduced.
900
c mxstep iwork(6) maximum number of (internally defined) steps
901
c allowed during one call to the solver.
902
c the default value is 500.
904
c mxhnil iwork(7) maximum number of messages printed (per problem)
905
c warning that t + h = t on a step (h = step size).
906
c this must be positive to result in a non-default
907
c value. the default value is 10.
908
c-----------------------------------------------------------------------
911
c as optional additional output from lsoibt, the variables listed
912
c below are quantities related to the performance of lsoibt
913
c which are available to the user. these are communicated by way of
914
c the work arrays, but also have internal mnemonic names as shown.
915
c except where stated otherwise, all of these outputs are defined
916
c on any successful return from lsoibt, and on any return with
917
c istate = -1, -2, -4, -5, -6, or -7. on a return with -3 (illegal
918
c input) or -8, they will be unchanged from their existing values
919
c (if any), except possibly for tolsf, lenrw, and leniw.
920
c on any error return, outputs relevant to the error will be defined,
923
c name location meaning
925
c hu rwork(11) the step size in t last used (successfully).
927
c hcur rwork(12) the step size to be attempted on the next step.
929
c tcur rwork(13) the current value of the independent variable
930
c which the solver has actually reached, i.e. the
931
c current internal mesh point in t. on output, tcur
932
c will always be at least as far as the argument
933
c t, but may be farther (if interpolation was done).
935
c tolsf rwork(14) a tolerance scale factor, greater than 1.0,
936
c computed when a request for too much accuracy was
937
c detected (istate = -3 if detected at the start of
938
c the problem, istate = -2 otherwise). if itol is
939
c left unaltered but rtol and atol are uniformly
940
c scaled up by a factor of tolsf for the next call,
941
c then the solver is deemed likely to succeed.
942
c (the user may also ignore tolsf and alter the
943
c tolerance parameters in any other way appropriate.)
945
c nst iwork(11) the number of steps taken for the problem so far.
947
c nre iwork(12) the number of residual evaluations (res calls)
948
c for the problem so far.
950
c nje iwork(13) the number of jacobian evaluations (each involving
951
c an evaluation of a and dr/dy) for the problem so
952
c far. this equals the number of calls to adda and
953
c (if miter = 1) to jac, and the number of matrix
954
c l-u decompositions.
956
c nqu iwork(14) the method order last used (successfully).
958
c nqcur iwork(15) the order to be attempted on the next step.
960
c imxer iwork(16) the index of the component of largest magnitude in
961
c the weighted local error vector ( e(i)/ewt(i) ),
962
c on an error return with istate = -4 or -5.
964
c lenrw iwork(17) the length of rwork actually required.
965
c this is defined on normal returns and on an illegal
966
c input return for insufficient storage.
968
c leniw iwork(18) the length of iwork actually required.
969
c this is defined on normal returns and on an illegal
970
c input return for insufficient storage.
973
c the following two arrays are segments of the rwork array which
974
c may also be of interest to the user as optional outputs.
975
c for each array, the table below gives its internal name,
976
c its base address in rwork, and its description.
978
c name base address description
980
c yh 21 the nordsieck history array, of size nyh by
981
c (nqcur + 1), where nyh is the initial value
982
c of neq. for j = 0,1,...,nqcur, column j+1
983
c of yh contains hcur**j/factorial(j) times
984
c the j-th derivative of the interpolating
985
c polynomial currently representing the solution,
986
c evaluated at t = tcur.
988
c acor lenrw-neq+1 array of size neq used for the accumulated
989
c corrections on each step, scaled on output to
990
c represent the estimated local error in y on the
991
c last step. this is the vector e in the descrip-
992
c tion of the error control. it is defined only
993
c on a return from lsoibt with istate = 2.
995
c-----------------------------------------------------------------------
996
c part ii. other routines callable.
998
c the following are optional calls which the user may make to
999
c gain additional capabilities in conjunction with lsoibt.
1000
c (the routines xsetun and xsetf are designed to conform to the
1001
c slatec error handling package.)
1003
c form of call function
1004
c call xsetun(lun) set the logical unit number, lun, for
1005
c output of messages from lsoibt, if
1006
c the default is not desired.
1007
c the default value of lun is 6.
1009
c call xsetf(mflag) set a flag to control the printing of
1010
c messages by lsoibt.
1011
c mflag = 0 means do not print. (danger..
1012
c this risks losing valuable information.)
1013
c mflag = 1 means print (the default).
1015
c either of the above calls may be made at
1016
c any time and will take effect immediately.
1018
c call srcom(rsav,isav,job) saves and restores the contents of
1019
c the internal common blocks used by
1020
c lsoibt (see part iii below).
1021
c rsav must be a real array of length 218
1022
c or more, and isav must be an integer
1023
c array of length 41 or more.
1024
c job=1 means save common into rsav/isav.
1025
c job=2 means restore common from rsav/isav.
1026
c srcom is useful if one is
1027
c interrupting a run and restarting
1028
c later, or alternating between two or
1029
c more problems solved with lsoibt.
1031
c call intdy(,,,,,) provide derivatives of y, of various
1032
c (see below) orders, at a specified point t, if
1033
c desired. it may be called only after
1034
c a successful return from lsoibt.
1036
c the detailed instructions for using intdy are as follows.
1037
c the form of the call is..
1039
c call intdy (t, k, rwork(21), nyh, dky, iflag)
1041
c the input parameters are..
1043
c t = value of independent variable where answers are desired
1044
c (normally the same as the t last returned by lsoibt).
1045
c for valid results, t must lie between tcur - hu and tcur.
1046
c (see optional outputs for tcur and hu.)
1047
c k = integer order of the derivative desired. k must satisfy
1048
c 0 .le. k .le. nqcur, where nqcur is the current order
1049
c (see optional outputs). the capability corresponding
1050
c to k = 0, i.e. computing y(t), is already provided
1051
c by lsoibt directly. since nqcur .ge. 1, the first
1052
c derivative dy/dt is always available with intdy.
1053
c rwork(21) = the base address of the history array yh.
1054
c nyh = column length of yh, equal to the initial value of neq.
1056
c the output parameters are..
1058
c dky = a real array of length neq containing the computed value
1059
c of the k-th derivative of y(t).
1060
c iflag = integer flag, returned as 0 if k and t were legal,
1061
c -1 if k was illegal, and -2 if t was illegal.
1062
c on an error return, a message is also written.
1063
c-----------------------------------------------------------------------
1064
c part iii. common blocks.
1066
c if lsoibt is to be used in an overlay situation, the user
1067
c must declare, in the primary overlay, the variables in..
1068
c (1) the call sequence to lsoibt,
1069
c (2) the two internal common blocks
1070
c /ls0001/ of length 257 (218 double precision words
1071
c followed by 39 integer words),
1072
c /eh0001/ of length 2 (integer words).
1074
c if lsoibt is used on a system in which the contents of internal
1075
c common blocks are not preserved between calls, the user should
1076
c declare the above two common blocks in his main program to insure
1077
c that their contents are preserved.
1079
c if the solution of a given problem by lsoibt is to be interrupted
1080
c and then later continued, such as when restarting an interrupted run
1081
c or alternating between two or more problems, the user should save,
1082
c following the return from the last lsoibt call prior to the
1083
c interruption, the contents of the call sequence variables and the
1084
c internal common blocks, and later restore these values before the
1085
c next lsoibt call for that problem. to save and restore the common
1086
c blocks, use subroutine srcom (see part ii above).
1088
c-----------------------------------------------------------------------
1089
c part iv. optionally replaceable solver routines.
1091
c below are descriptions of two routines in the lsoibt package which
1092
c relate to the measurement of errors. either routine can be
1093
c replaced by a user-supplied version, if desired. however, since such
1094
c a replacement may have a major impact on performance, it should be
1095
c done only when absolutely necessary, and only with great caution.
1096
c (note.. the means by which the package version of a routine is
1097
c superseded by the user-s version may be system-dependent.)
1100
c the following subroutine is called just before each internal
1101
c integration step, and sets the array of error weights, ewt, as
1102
c described under itol/rtol/atol above..
1103
c subroutine ewset (neq, itol, rtol, atol, ycur, ewt)
1104
c where neq, itol, rtol, and atol are as in the lsoibt call sequence,
1105
c ycur contains the current dependent variable vector, and
1106
c ewt is the array of weights set by ewset.
1108
c if the user supplies this subroutine, it must return in ewt(i)
1109
c (i = 1,...,neq) a positive quantity suitable for comparing errors
1110
c in y(i) to. the ewt array returned by ewset is passed to the
1111
c vnorm routine (see below), and also used by lsoibt in the computation
1112
c of the optional output imxer, the diagonal jacobian approximation,
1113
c and the increments for difference quotient jacobians.
1115
c in the user-supplied version of ewset, it may be desirable to use
1116
c the current values of derivatives of y. derivatives up to order nq
1117
c are available from the history array yh, described above under
1118
c optional outputs. in ewset, yh is identical to the ycur array,
1119
c extended to nq + 1 columns with a column length of nyh and scale
1120
c factors of h**j/factorial(j). on the first call for the problem,
1121
c given by nst = 0, nq is 1 and h is temporarily set to 1.0.
1122
c the quantities nq, nyh, h, and nst can be obtained by including
1123
c in ewset the statements..
1124
c double precision h, rls
1125
c common /ls0001/ rls(218),ils(39)
1130
c thus, for example, the current value of dy/dt can be obtained as
1131
c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is
1132
c unnecessary when nst = 0).
1135
c the following is a real function routine which computes the weighted
1136
c root-mean-square norm of a vector v..
1137
c d = vnorm (n, v, w)
1139
c n = the length of the vector,
1140
c v = real array of length n containing the vector,
1141
c w = real array of length n containing weights,
1142
c d = sqrt( (1/n) * sum(v(i)*w(i))**2 ).
1143
c vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where
1144
c ewt is as set by subroutine ewset.
1146
c if the user supplies this function, it should return a non-negative
1147
c value of vnorm suitable for use in the error control in lsoibt.
1148
c none of the arguments should be altered by vnorm.
1149
c for example, a user-supplied vnorm routine might..
1150
c -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
1151
c -ignore some components of v in the norm, with the effect of
1152
c suppressing the error control on those components of y.
1153
c-----------------------------------------------------------------------
1154
c-----------------------------------------------------------------------
1155
c other routines in the lsoibt package.
1157
c in addition to subroutine lsoibt, the lsoibt package includes the
1158
c following subroutines and function routines..
1159
c aigbt computes the initial value of the vector
1160
c dy/dt = inverse(a) * g
1161
c intdy computes an interpolated value of the y vector at t = tout.
1162
c stodi is the core integrator, which does one step of the
1163
c integration and the associated error control.
1164
c cfode sets all method coefficients and test constants.
1165
c ewset sets the error weight vector ewt before each step.
1166
c vnorm computes the weighted r.m.s. norm of a vector.
1167
c srcom is a user-callable routine to save and restore
1168
c the contents of the internal common blocks.
1169
c pjibt computes and preprocesses the jacobian matrix
1170
c and the newton iteration matrix p.
1171
c slsbt manages solution of linear system in chord iteration.
1172
c decbt and solbt are routines for solving block-tridiagonal
1173
c systems of linear algebraic equations.
1174
c dgefa and dgesl are routines from linpack for solving full
1175
c systems of linear algebraic equations.
1176
c daxpy, dscal, idamax, and ddot are basic linear algebra modules
1177
c (blas) used by the above linpack routines.
1178
c d1mach computes the unit roundoff in a machine-independent manner.
1179
c xerrwv, xsetun, and xsetf handle the printing of all error
1180
c messages and warnings. xerrwv is machine-dependent.
1181
c note.. vnorm, idamax, ddot, and d1mach are function routines.
1182
c all the others are subroutines.
1184
c the intrinsic and external routines used by lsoibt are.. dabs,
1185
c dmax1, dmin1, dfloat, iabs, max0, min0, mod, dsign, dsqrt, and write.
1187
c a block data subprogram is also included with the package,
1188
c for loading some of the variables in internal common.
1190
c-----------------------------------------------------------------------
1191
c the following card is for optimized compilation on llnl compilers.
1193
c-----------------------------------------------------------------------
1194
external pjibt, slsbt
1195
integer illin, init, lyh, lewt, lacor, lsavr, lwm, liwm,
1196
1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
1197
integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1198
1 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
1199
integer i, i1, i2, ier, iflag, imxer, ires, kgo,
1200
1 leniw, lenrw, lenwm, lp, lyd0, mb, mord, mxhnl0, mxstp0, nb
1201
double precision rowns,
1202
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1203
double precision atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
1204
1 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0,
1208
c-----------------------------------------------------------------------
1209
c the following internal common block contains
1210
c (a) variables which are local to any subroutine but whose values must
1211
c be preserved between calls to the routine (own variables), and
1212
c (b) variables which are communicated between subroutines.
1213
c block ls0001 is shared by the lsoibt, lsodi, and lsode packages.
1214
c the structure of ls0001 is as follows.. all real variables are
1215
c listed first, followed by all integers. within each type, the
1216
c variables are grouped with those local to subroutine lsoibt first,
1217
c then those local to subroutine stodi, and finally those used
1218
c for communication. the block is declared in subroutines
1219
c lsoibt, intdy, stodi, pjibt, and slsbt. groups of variables are
1220
c replaced by dummy arrays in the common declarations in routines
1221
c where those variables are not used.
1222
c-----------------------------------------------------------------------
1223
common /ls0001/ rowns(209),
1224
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1225
2 illin, init, lyh, lewt, lacor, lsavr, lwm, liwm,
1226
3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
1227
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1228
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
1230
data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1231
c-----------------------------------------------------------------------
1233
c this code block is executed on every call.
1234
c it tests istate and itask for legality and branches appropriately.
1235
c if istate .gt. 1 but the flag init shows that initialization has
1236
c not yet been done, an error return occurs.
1237
c if istate = 0 or 1 and tout = t, jump to block g and return
1239
c-----------------------------------------------------------------------
1240
if (istate .lt. 0 .or. istate .gt. 3) go to 601
1241
if (itask .lt. 1 .or. itask .gt. 5) go to 602
1242
if (istate .le. 1) go to 10
1243
if (init .eq. 0) go to 603
1244
if (istate .eq. 2) go to 200
1247
if (tout .eq. t) go to 430
1249
c-----------------------------------------------------------------------
1251
c the next code block is executed for the initial call (istate = 0 or 1)
1252
c or for a continuation call with parameter changes (istate = 3).
1253
c it contains checking of all inputs and various initializations.
1255
c first check legality of the non-optional inputs neq, itol, iopt,
1257
c-----------------------------------------------------------------------
1258
if (neq(1) .le. 0) go to 604
1259
if (istate .le. 1) go to 25
1260
if (neq(1) .gt. n) go to 605
1262
if (itol .lt. 1 .or. itol .gt. 4) go to 606
1263
if (iopt .lt. 0 .or. iopt .gt. 1) go to 607
1265
miter = mf - 10*meth
1266
if (meth .lt. 1 .or. meth .gt. 2) go to 608
1267
if (miter .lt. 1 .or. miter .gt. 2) go to 608
1270
if (mb .lt. 1 .or. mb .gt. n) go to 609
1271
if (nb .lt. 4) go to 610
1272
if (mb*nb .ne. n) go to 609
1273
c next process and check the optional inputs. --------------------------
1274
if (iopt .eq. 1) go to 40
1278
if (istate .le. 1) h0 = 0.0d0
1282
40 maxord = iwork(5)
1283
if (maxord .lt. 0) go to 611
1284
if (maxord .eq. 0) maxord = 100
1285
maxord = min0(maxord,mord(meth))
1287
if (mxstep .lt. 0) go to 612
1288
if (mxstep .eq. 0) mxstep = mxstp0
1290
if (mxhnil .lt. 0) go to 613
1291
if (mxhnil .eq. 0) mxhnil = mxhnl0
1292
if (istate .gt. 1) go to 50
1294
if ((tout - t)*h0 .lt. 0.0d0) go to 614
1296
if (hmax .lt. 0.0d0) go to 615
1298
if (hmax .gt. 0.0d0) hmxi = 1.0d0/hmax
1300
if (hmin .lt. 0.0d0) go to 616
1301
c-----------------------------------------------------------------------
1302
c set work array pointers and check lengths lrw and liw.
1303
c pointers to segments of rwork and iwork are named by prefixing l to
1304
c the name of the segment. e.g., the segment yh starts at rwork(lyh).
1305
c segments of rwork (in order) are denoted yh, wm, ewt, savr, acor.
1306
c-----------------------------------------------------------------------
1308
if (istate .le. 1) nyh = n
1309
lwm = lyh + (maxord + 1)*nyh
1310
lenwm = 3*mb*mb*nb + 2
1314
lenrw = lacor + n - 1
1319
if (lenrw .gt. lrw) go to 617
1320
if (leniw .gt. liw) go to 618
1321
c check rtol and atol for legality. ------------------------------------
1325
if (itol .ge. 3) rtoli = rtol(i)
1326
if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1327
if (rtoli .lt. 0.0d0) go to 619
1328
if (atoli .lt. 0.0d0) go to 620
1330
if (istate .le. 1) go to 100
1331
c if istate = 3, set flag to signal parameter changes to stodi. --------
1333
if (nq .le. maxord) go to 90
1334
c maxord was reduced below nq. copy yh(*,maxord+2) into ydoti.---------
1336
80 ydoti(i) = rwork(i+lwm-1)
1337
c reload wm(1) = rwork(lwm), since lwm may have changed. ---------------
1338
90 rwork(lwm) = dsqrt(uround)
1339
if (n .eq. nyh) go to 200
1340
c neq was reduced. zero part of yh to avoid undefined references. -----
1342
i2 = lyh + (maxord + 1)*nyh - 1
1343
if (i1 .gt. i2) go to 200
1347
c-----------------------------------------------------------------------
1349
c the next block is for the initial call only (istate = 0 or 1).
1350
c it contains all remaining initializations, the call to aigbt
1351
c (if istate = 1), and the calculation of the initial step size.
1352
c the error weights in ewt are inverted after being loaded.
1353
c-----------------------------------------------------------------------
1354
100 uround = d1mach(4)
1356
if (itask .ne. 4 .and. itask .ne. 5) go to 105
1358
if ((tcrit - tout)*(tout - t) .lt. 0.0d0) go to 625
1359
if (h0 .ne. 0.0d0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d0)
1362
rwork(lwm) = dsqrt(uround)
1374
c compute initial dy/dt, if necessary, and load it and initial y into yh
1377
if ( istate .eq. 1 ) go to 120
1378
c lsoibt must compute initial dy/dt (lyd0 points to yh(*,2)). ----------
1379
call aigbt( res, adda, neq, t, y, rwork(lyd0),
1380
1 mb, nb, rwork(lp), iwork(21), ier )
1382
if (ier) 560,110,565
1385
115 rwork(i+lyh-1) = y(i)
1387
c initial dy/dt has been supplied. -------------------------------------
1389
rwork(i+lyh-1) = y(i)
1390
125 rwork(i+lyd0-1) = ydoti(i)
1391
c load and invert the ewt array. (h is temporarily set to 1.0.) -------
1395
call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1397
if (rwork(i+lewt-1) .le. 0.0d0) go to 621
1398
135 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1399
c-----------------------------------------------------------------------
1400
c the coding below computes the step size, h0, to be attempted on the
1401
c first step, unless the user has supplied a value for this.
1402
c first check that tout - t differs significantly from zero.
1403
c a scalar tolerance quantity tol is computed, as max(rtol(i))
1404
c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted
1405
c so as to be between 100*uround and 1.0e-3.
1406
c then the computed value h0 is given by..
1408
c h0**2 = tol / ( w0**-2 + (1/neq) * sum ( ydot(i)/ywt(i) )**2 )
1410
c where w0 = max ( abs(t), abs(tout) ),
1411
c ydot(i) = i-th component of initial value of dy/dt,
1412
c ywt(i) = ewt(i)/tol (a weight for y(i)).
1413
c the sign of h0 is inferred from the initial values of tout and t.
1414
c-----------------------------------------------------------------------
1415
if (h0 .ne. 0.0d0) go to 180
1416
tdist = dabs(tout - t)
1417
w0 = dmax1(dabs(t),dabs(tout))
1418
if (tdist .lt. 2.0d0*uround*w0) go to 622
1420
if (itol .le. 2) go to 145
1422
140 tol = dmax1(tol,rtol(i))
1423
145 if (tol .gt. 0.0d0) go to 160
1426
if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1428
if (ayi .ne. 0.0d0) tol = dmax1(tol,atoli/ayi)
1430
160 tol = dmax1(tol,100.0d0*uround)
1431
tol = dmin1(tol,0.001d0)
1432
sum = vnorm (n, rwork(lyd0), rwork(lewt))
1433
sum = 1.0d0/(tol*w0*w0) + tol*sum**2
1434
h0 = 1.0d0/dsqrt(sum)
1435
h0 = dmin1(h0,tdist)
1436
h0 = dsign(h0,tout-t)
1437
c adjust h0 if necessary to meet hmax bound. ---------------------------
1438
180 rh = dabs(h0)*hmxi
1439
if (rh .gt. 1.0d0) h0 = h0/rh
1440
c load h with h0 and scale yh(*,2) by h0. ------------------------------
1443
190 rwork(i+lyd0-1) = h0*rwork(i+lyd0-1)
1445
c-----------------------------------------------------------------------
1447
c the next code block is for continuation calls only (istate = 2 or 3)
1448
c and is to check stop conditions before taking a step.
1449
c-----------------------------------------------------------------------
1451
go to (210, 250, 220, 230, 240), itask
1452
210 if ((tn - tout)*h .lt. 0.0d0) go to 250
1453
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1454
if (iflag .ne. 0) go to 627
1457
220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
1458
if ((tp - tout)*h .gt. 0.0d0) go to 623
1459
if ((tn - tout)*h .lt. 0.0d0) go to 250
1461
230 tcrit = rwork(1)
1462
if ((tn - tcrit)*h .gt. 0.0d0) go to 624
1463
if ((tcrit - tout)*h .lt. 0.0d0) go to 625
1464
if ((tn - tout)*h .lt. 0.0d0) go to 245
1465
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1466
if (iflag .ne. 0) go to 627
1469
240 tcrit = rwork(1)
1470
if ((tn - tcrit)*h .gt. 0.0d0) go to 624
1471
245 hmx = dabs(tn) + dabs(h)
1472
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1474
tnext = tn + h*(1.0d0 + 4.0d0*uround)
1475
if ((tnext - tcrit)*h .le. 0.0d0) go to 250
1476
h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1477
if (istate .eq. 2) jstart = -2
1478
c-----------------------------------------------------------------------
1480
c the next block is normally executed for all calls and contains
1481
c the call to the one-step core integrator stodi.
1483
c this is a looping point for the integration steps.
1485
c first check for too many steps being taken, update ewt (if not at
1486
c start of problem), check for too much accuracy being requested, and
1487
c check for h below the roundoff level in t.
1488
c-----------------------------------------------------------------------
1490
if ((nst-nslast) .ge. mxstep) go to 500
1491
call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1493
if (rwork(i+lewt-1) .le. 0.0d0) go to 510
1494
260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1495
270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt))
1496
if (tolsf .le. 1.0d0) go to 280
1498
if (nst .eq. 0) go to 626
1500
280 if ((tn + h) .ne. tn) go to 290
1502
if (nhnil .gt. mxhnil) go to 290
1503
call xerrwv(50hlsoibt-- warning..internal t (=r1) and h (=r2) are,
1504
1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1506
1 60h such that in the machine, t + h = t on the next step ,
1507
1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1508
call xerrwv(50h (h = step size). solver will continue anyway,
1509
1 50, 101, 0, 0, 0, 0, 2, tn, h)
1510
if (nhnil .lt. mxhnil) go to 290
1511
call xerrwv(50hlsoibt-- above warning has been issued i1 times. ,
1512
1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1513
call xerrwv(50h it will not be issued again for this problem,
1514
1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1516
c-----------------------------------------------------------------------
1517
c call stodi(neq,y,yh,nyh,yh1,ewt,savf,savr,acor,wm,iwm,res,
1518
c adda,jac,pjibt,slsbt)
1519
c note... savf in stodi occupies the same space as ydoti in lsoibt.
1520
c-----------------------------------------------------------------------
1521
call stodi (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1522
1 ydoti, rwork(lsavr), rwork(lacor), rwork(lwm),
1523
2 iwork(liwm), res, adda, jac, pjibt, slsbt )
1525
go to (300, 530, 540, 400, 550), kgo
1527
c kgo = 1,success. 2,error test failure. 3,convergence failure.
1528
c 4,res ordered return. 5,res returned error.
1529
c-----------------------------------------------------------------------
1531
c the following block handles the case of a successful return from the
1532
c core integrator (kflag = 0). test for stop conditions.
1533
c-----------------------------------------------------------------------
1535
go to (310, 400, 330, 340, 350), itask
1536
c itask = 1. if tout has been reached, interpolate. -------------------
1537
310 if ((tn - tout)*h .lt. 0.0d0) go to 250
1538
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1541
c itask = 3. jump to exit if tout was reached. ------------------------
1542
330 if ((tn - tout)*h .ge. 0.0d0) go to 400
1544
c itask = 4. see if tout or tcrit was reached. adjust h if necessary.
1545
340 if ((tn - tout)*h .lt. 0.0d0) go to 345
1546
call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1549
345 hmx = dabs(tn) + dabs(h)
1550
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1552
tnext = tn + h*(1.0d0 + 4.0d0*uround)
1553
if ((tnext - tcrit)*h .le. 0.0d0) go to 250
1554
h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1557
c itask = 5. see if tcrit was reached and jump to exit. ---------------
1558
350 hmx = dabs(tn) + dabs(h)
1559
ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1560
c-----------------------------------------------------------------------
1562
c the following block handles all successful returns from lsoibt.
1563
c if itask .ne. 1, y is loaded from yh and t is set accordingly.
1564
c istate is set to 2, the illegal input counter is zeroed, and the
1565
c optional outputs are loaded into the work arrays before returning. if
1566
c istate = 0 or 1 and tout = t, there is a return with no action taken,
1567
c except that if this has happened repeatedly, the run is terminated.
1568
c-----------------------------------------------------------------------
1570
410 y(i) = rwork(i+lyh-1)
1572
if (itask .ne. 4 .and. itask .ne. 5) go to 420
1575
if ( kflag .eq. -3 ) istate = 3
1587
430 ntrep = ntrep + 1
1588
if (ntrep .lt. 5) return
1590
1 60hlsoibt-- repeated calls with istate= 0 or 1 and tout= t(=r1),
1591
1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0)
1593
c-----------------------------------------------------------------------
1595
c the following block handles all unsuccessful returns other than
1596
c those for illegal input. first the error message routine is called.
1597
c if there was an error test or convergence test failure, imxer is set.
1598
c then y is loaded from yh, t is set to tn, and the illegal input
1599
c counter illin is set to 0. the optional outputs are loaded into
1600
c the work arrays before returning.
1601
c-----------------------------------------------------------------------
1602
c the maximum number of steps was taken before reaching tout. ----------
1603
500 call xerrwv(50hlsoibt-- at current t (=r1), mxstep (=i1) steps ,
1604
1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1605
call xerrwv(50h taken on this call before reaching tout ,
1606
1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
1609
c ewt(i) .le. 0.0 for some i (not at start of problem). ----------------
1610
510 ewti = rwork(lewt+i-1)
1611
call xerrwv(50hlsoibt-- at t (=r1), ewt(i1) has become r2 .le. 0.,
1612
1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1615
c too much accuracy requested for machine precision. -------------------
1616
520 call xerrwv(50hlsoibt-- at t (=r1), too much accuracy requested ,
1617
1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1618
call xerrwv(50h for precision of machine.. see tolsf (=r2) ,
1619
1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1623
c kflag = -1. error test failed repeatedly or with abs(h) = hmin. -----
1624
530 call xerrwv(50hlsoibt-- at t(=r1) and step size h(=r2), the error,
1625
1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1626
call xerrwv(50h test failed repeatedly or with abs(h) = hmin,
1627
1 50, 204, 0, 0, 0, 0, 2, tn, h)
1630
c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ----
1631
540 call xerrwv(50hlsoibt-- at t (=r1) and step size h (=r2), the ,
1632
1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1633
call xerrwv(50h corrector convergence failed repeatedly ,
1634
1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1635
call xerrwv(30h or with abs(h) = hmin ,
1636
1 30, 205, 0, 0, 0, 0, 2, tn, h)
1639
c ires = 3 returned by res, despite retries by stodi. ------------------
1640
550 call xerrwv(50hlsoibt-- at t (=r1) residual routine returned ,
1641
1 50, 206, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1642
call xerrwv(40h error ires = 3 repeatedly ,
1643
1 40, 206, 0, 0, 0, 0, 1, tn, 0.0d0)
1646
c aigbt failed because a diagonal block of a-matrix was singular. ------
1649
1 60hlsoibt-- attempt to initialize dy/dt failed.. matrix a has a,
1650
1 60, 207, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1651
call xerrwv(50h singular diagonal block, block no. = (i1) ,
1652
2 50, 207, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
1655
c aigbt failed because res set ires to 2 or 3. -------------------------
1656
565 call xerrwv(50hlsoibt-- attempt to initialize dy/dt failed ,
1657
1 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1658
call xerrwv(50h because residual routine set its error flag ,
1659
1 50, 208, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1660
call xerrwv(20h to ires = (i1),
1661
1 20, 208, 0, 1, ier, 0, 0, 0.0d0, 0.0d0)
1664
c compute imxer if relevant. -------------------------------------------
1668
size = dabs(rwork(i+lacor-1)*rwork(i+lewt-1))
1669
if (big .ge. size) go to 575
1674
c compute residual if relevant. ----------------------------------------
1675
580 lyd0 = lyh + nyh
1677
rwork(i+lsavr-1) = rwork(i+lyd0-1)/h
1678
585 y(i) = rwork(i+lyh-1)
1680
call res ( neq, tn, y, rwork(lsavr), ydoti, ires )
1682
if ( ires .le. 1 ) go to 595
1683
call xerrwv(50hlsoibt-- residual routine set its flag ires ,
1684
1 50, 210, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1685
call xerrwv(50h to (i1) when called for final output. ,
1686
1 50, 210, 0, 1, ires, 0, 0, 0.0d0, 0.0d0)
1688
c set y vector, t, illin, and optional outputs. ------------------------
1690
592 y(i) = rwork(i+lyh-1)
1702
c-----------------------------------------------------------------------
1704
c the following block handles all error returns due to illegal input
1705
c (istate = -3), as detected before calling the core integrator.
1706
c first the error message routine is called. then if there have been
1707
c 5 consecutive such returns just before this call to the solver,
1708
c the run is halted.
1709
c-----------------------------------------------------------------------
1710
601 call xerrwv(30hlsoibt-- istate (=i1) illegal ,
1711
1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
1713
602 call xerrwv(30hlsoibt-- itask (=i1) illegal ,
1714
1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
1716
603 call xerrwv(50hlsoibt-- istate .gt. 1 but lsoibt not initialized ,
1717
1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1719
604 call xerrwv(30hlsoibt-- neq (=i1) .lt. 1 ,
1720
1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
1722
605 call xerrwv(50hlsoibt-- istate = 3 and neq increased (i1 to i2) ,
1723
1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
1725
606 call xerrwv(30hlsoibt-- itol (=i1) illegal ,
1726
1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
1728
607 call xerrwv(30hlsoibt-- iopt (=i1) illegal ,
1729
1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
1731
608 call xerrwv(30hlsoibt-- mf (=i1) illegal ,
1732
1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
1734
609 call xerrwv(40hlsoibt-- mb (=i1) or nb (=i2) illegal ,
1735
1 50, 9, 0, 2, mb, nb, 0, 0.0d0, 0.0d0)
1737
610 call xerrwv(40hlsoibt-- nb(=i1) illegal.. .lt. 4 ,
1738
1 50, 10, 0, 1, nb, 0, 0, 0.0d0, 0.0d0)
1740
611 call xerrwv(30hlsoibt-- maxord (=i1) .lt. 0 ,
1741
1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
1743
612 call xerrwv(30hlsoibt-- mxstep (=i1) .lt. 0 ,
1744
1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
1746
613 call xerrwv(30hlsoibt-- mxhnil (=i1) .lt. 0 ,
1747
1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1749
614 call xerrwv(40hlsoibt-- tout (=r1) behind t (=r2) ,
1750
1 40, 14, 0, 0, 0, 0, 2, tout, t)
1751
call xerrwv(50h integration direction is given by h0 (=r1) ,
1752
1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
1754
615 call xerrwv(30hlsoibt-- hmax (=r1) .lt. 0.0 ,
1755
1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
1757
616 call xerrwv(30hlsoibt-- hmin (=r1) .lt. 0.0 ,
1758
1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
1761
1 60hlsoibt-- rwork length needed, lenrw (=i1), exceeds lrw (=i2),
1762
1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
1765
1 60hlsoibt-- iwork length needed, leniw (=i1), exceeds liw (=i2),
1766
1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
1768
619 call xerrwv(40hlsoibt-- rtol(=i1) is r1 .lt. 0.0 ,
1769
1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
1771
620 call xerrwv(40hlsoibt-- atol(=i1) is r1 .lt. 0.0 ,
1772
1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
1774
621 ewti = rwork(lewt+i-1)
1775
call xerrwv(40hlsoibt-- ewt(=i1) is r1 .le. 0.0 ,
1776
1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
1779
1 60hlsoibt-- tout (=r1) too close to t(=r2) to start integration,
1780
1 60, 22, 0, 0, 0, 0, 2, tout, t)
1783
1 60hlsoibt-- itask = i1 and tout (=r1) behind tcur - hu (= r2) ,
1784
1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1787
1 60hlsoibt-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ,
1788
1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1791
1 60hlsoibt-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ,
1792
1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1794
626 call xerrwv(50hlsoibt-- at start of problem, too much accuracy ,
1795
1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1797
1 60h requested for precision of machine.. see tolsf (=r1) ,
1798
1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
1801
627 call xerrwv(50hlsoibt-- trouble from intdy. itask = i1, tout = r1,
1802
1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
1804
700 if (illin .eq. 5) go to 710
1808
710 call xerrwv(50hlsoibt-- repeated occurrences of illegal input ,
1809
1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1811
800 call xerrwv(50hlsoibt-- run aborted.. apparent infinite loop ,
1812
1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
1814
c----------------------- end of subroutine lsoibt ----------------------
added
removed
Lib/integrate/odepack/lsoibt.f
