Parent Directory | Revision Log

Revision **1507** -
(**show annotations**)
(**download**)

*Wed Jun 27 11:25:37 2007 UTC*
(16 years, 11 months ago)
by *jpye*

File size: 118049 byte(s)

File size: 118049 byte(s)

Moving integrators to own directory, about to make them self-contained shared libraries.

1 | subroutine lsode (f, neq, y, t, tout, itol, rtol, atol, itask, |

2 | 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf) |

3 | external f, jac, xascwv |

4 | integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf |

5 | double precision y, t, tout, rtol, atol, rwork |

6 | dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw) |

7 | c----------------------------------------------------------------------- |

8 | c this is the march 30, 1987 version of |

9 | c lsode.. livermore solver for ordinary differential equations. |

10 | c this version is in double precision. |

11 | c |

12 | c lsode solves the initial value problem for stiff or nonstiff |

13 | c systems of first order ode-s, |

14 | c dy/dt = f(t,y) , or, in component form, |

15 | c dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq). |

16 | c lsode is a package based on the gear and gearb packages, and on the |

17 | c october 23, 1978 version of the tentative odepack user interface |

18 | c standard, with minor modifications. |

19 | c----------------------------------------------------------------------- |

20 | c reference.. |

21 | c alan c. hindmarsh, odepack, a systematized collection of ode |

22 | c solvers, in scientific computing, r. s. stepleman et al. (eds.), |

23 | c north-holland, amsterdam, 1983, pp. 55-64. |

24 | c----------------------------------------------------------------------- |

25 | c author and contact.. alan c. hindmarsh, |

26 | c computing and mathematics research div., l-316 |

27 | c lawrence livermore national laboratory |

28 | c livermore, ca 94550. |

29 | c----------------------------------------------------------------------- |

30 | c summary of usage. |

31 | c |

32 | c communication between the user and the lsode package, for normal |

33 | c situations, is summarized here. this summary describes only a subset |

34 | c of the full set of options available. see the full description for |

35 | c details, including optional communication, nonstandard options, |

36 | c and instructions for special situations. see also the example |

37 | c problem (with program and output) following this summary. |

38 | c |

39 | c a. first provide a subroutine of the form.. |

40 | c subroutine f (neq, t, y, ydot) |

41 | c dimension y(neq), ydot(neq) |

42 | c which supplies the vector function f by loading ydot(i) with f(i). |

43 | c |

44 | c b. next determine (or guess) whether or not the problem is stiff. |

45 | c stiffness occurs when the jacobian matrix df/dy has an eigenvalue |

46 | c whose real part is negative and large in magnitude, compared to the |

47 | c reciprocal of the t span of interest. if the problem is nonstiff, |

48 | c use a method flag mf = 10. if it is stiff, there are four standard |

49 | c choices for mf, and lsode requires the jacobian matrix in some form. |

50 | c this matrix is regarded either as full (mf = 21 or 22), |

51 | c or banded (mf = 24 or 25). in the banded case, lsode requires two |

52 | c half-bandwidth parameters ml and mu. these are, respectively, the |

53 | c widths of the lower and upper parts of the band, excluding the main |

54 | c diagonal. thus the band consists of the locations (i,j) with |

55 | c i-ml .le. j .le. i+mu, and the full bandwidth is ml+mu+1. |

56 | c |

57 | c c. if the problem is stiff, you are encouraged to supply the jacobian |

58 | c directly (mf = 21 or 24), but if this is not feasible, lsode will |

59 | c compute it internally by difference quotients (mf = 22 or 25). |

60 | c if you are supplying the jacobian, provide a subroutine of the form.. |

61 | c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) |

62 | c dimension y(neq), pd(nrowpd,neq) |

63 | c which supplies df/dy by loading pd as follows.. |

64 | c for a full jacobian (mf = 21), load pd(i,j) with df(i)/dy(j), |

65 | c the partial derivative of f(i) with respect to y(j). (ignore the |

66 | c ml and mu arguments in this case.) |

67 | c for a banded jacobian (mf = 24), load pd(i-j+mu+1,j) with |

68 | c df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of |

69 | c pd from the top down. |

70 | c in either case, only nonzero elements need be loaded. |

71 | c |

72 | c d. write a main program which calls subroutine lsode once for |

73 | c each point at which answers are desired. this should also provide |

74 | c for possible use of logical unit 6 for output of error messages |

75 | c by lsode. on the first call to lsode, supply arguments as follows.. |

76 | c f = name of subroutine for right-hand side vector f. |

77 | c this name must be declared external in calling program. |

78 | c neq = number of first order ode-s. |

79 | c y = array of initial values, of length neq. |

80 | c t = the initial value of the independent variable. |

81 | c tout = first point where output is desired (.ne. t). |

82 | c itol = 1 or 2 according as atol (below) is a scalar or array. |

83 | c rtol = relative tolerance parameter (scalar). |

84 | c atol = absolute tolerance parameter (scalar or array). |

85 | c the estimated local error in y(i) will be controlled so as |

86 | c to be roughly less (in magnitude) than |

87 | c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or |

88 | c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2. |

89 | c thus the local error test passes if, in each component, |

90 | c either the absolute error is less than atol (or atol(i)), |

91 | c or the relative error is less than rtol. |

92 | c use rtol = 0.0 for pure absolute error control, and |

93 | c use atol = 0.0 (or atol(i) = 0.0) for pure relative error |

94 | c control. caution.. actual (global) errors may exceed these |

95 | c local tolerances, so choose them conservatively. |

96 | c itask = 1 for normal computation of output values of y at t = tout. |

97 | c istate = integer flag (input and output). set istate = 1. |

98 | c iopt = 0 to indicate no optional inputs used. |

99 | c rwork = real work array of length at least.. |

100 | c 20 + 16*neq for mf = 10, |

101 | c 22 + 9*neq + neq**2 for mf = 21 or 22, |

102 | c 22 + 10*neq + (2*ml + mu)*neq for mf = 24 or 25. |

103 | c lrw = declared length of rwork (in user-s dimension). |

104 | c iwork = integer work array of length at least.. |

105 | c 20 for mf = 10, |

106 | c 20 + neq for mf = 21, 22, 24, or 25. |

107 | c if mf = 24 or 25, input in iwork(1),iwork(2) the lower |

108 | c and upper half-bandwidths ml,mu. |

109 | c liw = declared length of iwork (in user-s dimension). |

110 | c jac = name of subroutine for jacobian matrix (mf = 21 or 24). |

111 | c if used, this name must be declared external in calling |

112 | c program. if not used, pass a dummy name. |

113 | c mf = method flag. standard values are.. |

114 | c 10 for nonstiff (adams) method, no jacobian used. |

115 | c 21 for stiff (bdf) method, user-supplied full jacobian. |

116 | c 22 for stiff method, internally generated full jacobian. |

117 | c 24 for stiff method, user-supplied banded jacobian. |

118 | c 25 for stiff method, internally generated banded jacobian. |

119 | c note that the main program must declare arrays y, rwork, iwork, |

120 | c and possibly atol. |

121 | c |

122 | c e. the output from the first call (or any call) is.. |

123 | c y = array of computed values of y(t) vector. |

124 | c t = corresponding value of independent variable (normally tout). |

125 | c istate = 2 if lsode was successful, negative otherwise. |

126 | c -1 means excess work done on this call (perhaps wrong mf). |

127 | c -2 means excess accuracy requested (tolerances too small). |

128 | c -3 means illegal input detected (see printed message). |

129 | c -4 means repeated error test failures (check all inputs). |

130 | c -5 means repeated convergence failures (perhaps bad jacobian |

131 | c supplied or wrong choice of mf or tolerances). |

132 | c -6 means error weight became zero during problem. (solution |

133 | c component i vanished, and atol or atol(i) = 0.) |

134 | c |

135 | c f. to continue the integration after a successful return, simply |

136 | c reset tout and call lsode again. no other parameters need be reset. |

137 | c |

138 | c----------------------------------------------------------------------- |

139 | c example problem. |

140 | c |

141 | c the following is a simple example problem, with the coding |

142 | c needed for its solution by lsode. the problem is from chemical |

143 | c kinetics, and consists of the following three rate equations.. |

144 | c dy1/dt = -.04*y1 + 1.e4*y2*y3 |

145 | c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2 |

146 | c dy3/dt = 3.e7*y2**2 |

147 | c on the interval from t = 0.0 to t = 4.e10, with initial conditions |

148 | c y1 = 1.0, y2 = y3 = 0. the problem is stiff. |

149 | c |

150 | c the following coding solves this problem with lsode, using mf = 21 |

151 | c and printing results at t = .4, 4., ..., 4.e10. it uses |

152 | c itol = 2 and atol much smaller for y2 than y1 or y3 because |

153 | c y2 has much smaller values. |

154 | c at the end of the run, statistical quantities of interest are |

155 | c printed (see optional outputs in the full description below). |

156 | c |

157 | c external fex, jex |

158 | c double precision atol, rtol, rwork, t, tout, y |

159 | c dimension y(3), atol(3), rwork(58), iwork(23) |

160 | c neq = 3 |

161 | c y(1) = 1.d0 |

162 | c y(2) = 0.d0 |

163 | c y(3) = 0.d0 |

164 | c t = 0.d0 |

165 | c tout = .4d0 |

166 | c itol = 2 |

167 | c rtol = 1.d-4 |

168 | c atol(1) = 1.d-6 |

169 | c atol(2) = 1.d-10 |

170 | c atol(3) = 1.d-6 |

171 | c itask = 1 |

172 | c istate = 1 |

173 | c iopt = 0 |

174 | c lrw = 58 |

175 | c liw = 23 |

176 | c mf = 21 |

177 | c do 40 iout = 1,12 |

178 | c call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate, |

179 | c 1 iopt,rwork,lrw,iwork,liw,jex,mf) |

180 | c write(6,20)t,y(1),y(2),y(3) |

181 | c 20 format(7h at t =,e12.4,6h y =,3e14.6) |

182 | c if (istate .lt. 0) go to 80 |

183 | c 40 tout = tout*10.d0 |

184 | c write(6,60)iwork(11),iwork(12),iwork(13) |

185 | c 60 format(/12h no. steps =,i4,11h no. f-s =,i4,11h no. j-s =,i4) |

186 | c stop |

187 | c 80 write(6,90)istate |

188 | c 90 format(///22h error halt.. istate =,i3) |

189 | c stop |

190 | c end |

191 | c |

192 | c subroutine fex (neq, t, y, ydot) |

193 | c double precision t, y, ydot |

194 | c dimension y(3), ydot(3) |

195 | c ydot(1) = -.04d0*y(1) + 1.d4*y(2)*y(3) |

196 | c ydot(3) = 3.d7*y(2)*y(2) |

197 | c ydot(2) = -ydot(1) - ydot(3) |

198 | c return |

199 | c end |

200 | c |

201 | c subroutine jex (neq, t, y, ml, mu, pd, nrpd) |

202 | c double precision pd, t, y |

203 | c dimension y(3), pd(nrpd,3) |

204 | c pd(1,1) = -.04d0 |

205 | c pd(1,2) = 1.d4*y(3) |

206 | c pd(1,3) = 1.d4*y(2) |

207 | c pd(2,1) = .04d0 |

208 | c pd(2,3) = -pd(1,3) |

209 | c pd(3,2) = 6.d7*y(2) |

210 | c pd(2,2) = -pd(1,2) - pd(3,2) |

211 | c return |

212 | c end |

213 | c |

214 | c the output of this program (on a cdc-7600 in single precision) |

215 | c is as follows.. |

216 | c |

217 | c at t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02 |

218 | c at t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02 |

219 | c at t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01 |

220 | c at t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01 |

221 | c at t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01 |

222 | c at t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01 |

223 | c at t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01 |

224 | c at t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01 |

225 | c at t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01 |

226 | c at t = 4.0000e+08 y = 5.494529e-06 2.197824e-11 9.999945e-01 |

227 | c at t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01 |

228 | c at t = 4.0000e+10 y = -7.170586e-08 -2.868234e-13 1.000000e+00 |

229 | c |

230 | c no. steps = 330 no. f-s = 405 no. j-s = 69 |

231 | c----------------------------------------------------------------------- |

232 | c full description of user interface to lsode. |

233 | c |

234 | c the user interface to lsode consists of the following parts. |

235 | c |

236 | c i. the call sequence to subroutine lsode, which is a driver |

237 | c routine for the solver. this includes descriptions of both |

238 | c the call sequence arguments and of user-supplied routines. |

239 | c following these descriptions is a description of |

240 | c optional inputs available through the call sequence, and then |

241 | c a description of optional outputs (in the work arrays). |

242 | c |

243 | c ii. descriptions of other routines in the lsode package that may be |

244 | c (optionally) called by the user. these provide the ability to |

245 | c alter error message handling, save and restore the internal |

246 | c common, and obtain specified derivatives of the solution y(t). |

247 | c |

248 | c iii. descriptions of common blocks to be declared in overlay |

249 | c or similar environments, or to be saved when doing an interrupt |

250 | c of the problem and continued solution later. |

251 | c |

252 | c iv. description of two routines in the lsode package, either of |

253 | c which the user may replace with his own version, if desired. |

254 | c these relate to the measurement of errors. |

255 | c |

256 | c----------------------------------------------------------------------- |

257 | c part i. call sequence. |

258 | c |

259 | c the call sequence parameters used for input only are |

260 | c f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf, |

261 | c and those used for both input and output are |

262 | c y, t, istate. |

263 | c the work arrays rwork and iwork are also used for conditional and |

264 | c optional inputs and optional outputs. (the term output here refers |

265 | c to the return from subroutine lsode to the user-s calling program.) |

266 | c |

267 | c the legality of input parameters will be thoroughly checked on the |

268 | c initial call for the problem, but not checked thereafter unless a |

269 | c change in input parameters is flagged by istate = 3 on input. |

270 | c |

271 | c the descriptions of the call arguments are as follows. |

272 | c |

273 | c f = the name of the user-supplied subroutine defining the |

274 | c ode system. the system must be put in the first-order |

275 | c form dy/dt = f(t,y), where f is a vector-valued function |

276 | c of the scalar t and the vector y. subroutine f is to |

277 | c compute the function f. it is to have the form |

278 | c subroutine f (neq, t, y, ydot) |

279 | c dimension y(1), ydot(1) |

280 | c where neq, t, and y are input, and the array ydot = f(t,y) |

281 | c is output. y and ydot are arrays of length neq. |

282 | c (in the dimension statement above, 1 is a dummy |

283 | c dimension.. it can be replaced by any value.) |

284 | c subroutine f should not alter y(1),...,y(neq). |

285 | c f must be declared external in the calling program. |

286 | c |

287 | c subroutine f may access user-defined quantities in |

288 | c neq(2),... and/or in y(neq(1)+1),... if neq is an array |

289 | c (dimensioned in f) and/or y has length exceeding neq(1). |

290 | c see the descriptions of neq and y below. |

291 | c |

292 | c if quantities computed in the f routine are needed |

293 | c externally to lsode, an extra call to f should be made |

294 | c for this purpose, for consistent and accurate results. |

295 | c if only the derivative dy/dt is needed, use intdy instead. |

296 | c |

297 | c neq = the size of the ode system (number of first order |

298 | c ordinary differential equations). used only for input. |

299 | c neq may be decreased, but not increased, during the problem. |

300 | c if neq is decreased (with istate = 3 on input), the |

301 | c remaining components of y should be left undisturbed, if |

302 | c these are to be accessed in f and/or jac. |

303 | c |

304 | c normally, neq is a scalar, and it is generally referred to |

305 | c as a scalar in this user interface description. however, |

306 | c neq may be an array, with neq(1) set to the system size. |

307 | c (the lsode package accesses only neq(1).) in either case, |

308 | c this parameter is passed as the neq argument in all calls |

309 | c to f and jac. hence, if it is an array, locations |

310 | c neq(2),... may be used to store other integer data and pass |

311 | c it to f and/or jac. subroutines f and/or jac must include |

312 | c neq in a dimension statement in that case. |

313 | c |

314 | c y = a real array for the vector of dependent variables, of |

315 | c length neq or more. used for both input and output on the |

316 | c first call (istate = 1), and only for output on other calls. |

317 | c on the first call, y must contain the vector of initial |

318 | c values. on output, y contains the computed solution vector, |

319 | c evaluated at t. if desired, the y array may be used |

320 | c for other purposes between calls to the solver. |

321 | c |

322 | c this array is passed as the y argument in all calls to |

323 | c f and jac. hence its length may exceed neq, and locations |

324 | c y(neq+1),... may be used to store other real data and |

325 | c pass it to f and/or jac. (the lsode package accesses only |

326 | c y(1),...,y(neq).) |

327 | c |

328 | c t = the independent variable. on input, t is used only on the |

329 | c first call, as the initial point of the integration. |

330 | c on output, after each call, t is the value at which a |

331 | c computed solution y is evaluated (usually the same as tout). |

332 | c on an error return, t is the farthest point reached. |

333 | c |

334 | c tout = the next value of t at which a computed solution is desired. |

335 | c used only for input. |

336 | c |

337 | c when starting the problem (istate = 1), tout may be equal |

338 | c to t for one call, then should .ne. t for the next call. |

339 | c for the initial t, an input value of tout .ne. t is used |

340 | c in order to determine the direction of the integration |

341 | c (i.e. the algebraic sign of the step sizes) and the rough |

342 | c scale of the problem. integration in either direction |

343 | c (forward or backward in t) is permitted. |

344 | c |

345 | c if itask = 2 or 5 (one-step modes), tout is ignored after |

346 | c the first call (i.e. the first call with tout .ne. t). |

347 | c otherwise, tout is required on every call. |

348 | c |

349 | c if itask = 1, 3, or 4, the values of tout need not be |

350 | c monotone, but a value of tout which backs up is limited |

351 | c to the current internal t interval, whose endpoints are |

352 | c tcur - hu and tcur (see optional outputs, below, for |

353 | c tcur and hu). |

354 | c |

355 | c itol = an indicator for the type of error control. see |

356 | c description below under atol. used only for input. |

357 | c |

358 | c rtol = a relative error tolerance parameter, either a scalar or |

359 | c an array of length neq. see description below under atol. |

360 | c input only. |

361 | c |

362 | c atol = an absolute error tolerance parameter, either a scalar or |

363 | c an array of length neq. input only. |

364 | c |

365 | c the input parameters itol, rtol, and atol determine |

366 | c the error control performed by the solver. the solver will |

367 | c control the vector e = (e(i)) of estimated local errors |

368 | c in y, according to an inequality of the form |

369 | c rms-norm of ( e(i)/ewt(i) ) .le. 1, |

370 | c where ewt(i) = rtol(i)*abs(y(i)) + atol(i), |

371 | c and the rms-norm (root-mean-square norm) here is |

372 | c rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i)) |

373 | c is a vector of weights which must always be positive, and |

374 | c the values of rtol and atol should all be non-negative. |

375 | c the following table gives the types (scalar/array) of |

376 | c rtol and atol, and the corresponding form of ewt(i). |

377 | c |

378 | c itol rtol atol ewt(i) |

379 | c 1 scalar scalar rtol*abs(y(i)) + atol |

380 | c 2 scalar array rtol*abs(y(i)) + atol(i) |

381 | c 3 array scalar rtol(i)*abs(y(i)) + atol |

382 | c 4 array array rtol(i)*abs(y(i)) + atol(i) |

383 | c |

384 | c when either of these parameters is a scalar, it need not |

385 | c be dimensioned in the user-s calling program. |

386 | c |

387 | c if none of the above choices (with itol, rtol, and atol |

388 | c fixed throughout the problem) is suitable, more general |

389 | c error controls can be obtained by substituting |

390 | c user-supplied routines for the setting of ewt and/or for |

391 | c the norm calculation. see part iv below. |

392 | c |

393 | c if global errors are to be estimated by making a repeated |

394 | c run on the same problem with smaller tolerances, then all |

395 | c components of rtol and atol (i.e. of ewt) should be scaled |

396 | c down uniformly. |

397 | c |

398 | c itask = an index specifying the task to be performed. |

399 | c input only. itask has the following values and meanings. |

400 | c 1 means normal computation of output values of y(t) at |

401 | c t = tout (by overshooting and interpolating). |

402 | c 2 means take one step only and return. |

403 | c 3 means stop at the first internal mesh point at or |

404 | c beyond t = tout and return. |

405 | c 4 means normal computation of output values of y(t) at |

406 | c t = tout but without overshooting t = tcrit. |

407 | c tcrit must be input as rwork(1). tcrit may be equal to |

408 | c or beyond tout, but not behind it in the direction of |

409 | c integration. this option is useful if the problem |

410 | c has a singularity at or beyond t = tcrit. |

411 | c 5 means take one step, without passing tcrit, and return. |

412 | c tcrit must be input as rwork(1). |

413 | c |

414 | c note.. if itask = 4 or 5 and the solver reaches tcrit |

415 | c (within roundoff), it will return t = tcrit (exactly) to |

416 | c indicate this (unless itask = 4 and tout comes before tcrit, |

417 | c in which case answers at t = tout are returned first). |

418 | c |

419 | c istate = an index used for input and output to specify the |

420 | c the state of the calculation. |

421 | c |

422 | c on input, the values of istate are as follows. |

423 | c 1 means this is the first call for the problem |

424 | c (initializations will be done). see note below. |

425 | c 2 means this is not the first call, and the calculation |

426 | c is to continue normally, with no change in any input |

427 | c parameters except possibly tout and itask. |

428 | c (if itol, rtol, and/or atol are changed between calls |

429 | c with istate = 2, the new values will be used but not |

430 | c tested for legality.) |

431 | c 3 means this is not the first call, and the |

432 | c calculation is to continue normally, but with |

433 | c a change in input parameters other than |

434 | c tout and itask. changes are allowed in |

435 | c neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu, |

436 | c and any of the optional inputs except h0. |

437 | c (see iwork description for ml and mu.) |

438 | c note.. a preliminary call with tout = t is not counted |

439 | c as a first call here, as no initialization or checking of |

440 | c input is done. (such a call is sometimes useful for the |

441 | c purpose of outputting the initial conditions.) |

442 | c thus the first call for which tout .ne. t requires |

443 | c istate = 1 on input. |

444 | c |

445 | c on output, istate has the following values and meanings. |

446 | c 1 means nothing was done, as tout was equal to t with |

447 | c istate = 1 on input. (however, an internal counter was |

448 | c set to detect and prevent repeated calls of this type.) |

449 | c 2 means the integration was performed successfully. |

450 | c -1 means an excessive amount of work (more than mxstep |

451 | c steps) was done on this call, before completing the |

452 | c requested task, but the integration was otherwise |

453 | c successful as far as t. (mxstep is an optional input |

454 | c and is normally 500.) to continue, the user may |

455 | c simply reset istate to a value .gt. 1 and call again |

456 | c (the excess work step counter will be reset to 0). |

457 | c in addition, the user may increase mxstep to avoid |

458 | c this error return (see below on optional inputs). |

459 | c -2 means too much accuracy was requested for the precision |

460 | c of the machine being used. this was detected before |

461 | c completing the requested task, but the integration |

462 | c was successful as far as t. to continue, the tolerance |

463 | c parameters must be reset, and istate must be set |

464 | c to 3. the optional output tolsf may be used for this |

465 | c purpose. (note.. if this condition is detected before |

466 | c taking any steps, then an illegal input return |

467 | c (istate = -3) occurs instead.) |

468 | c -3 means illegal input was detected, before taking any |

469 | c integration steps. see written message for details. |

470 | c note.. if the solver detects an infinite loop of calls |

471 | c to the solver with illegal input, it will cause |

472 | c the run to stop. |

473 | c -4 means there were repeated error test failures on |

474 | c one attempted step, before completing the requested |

475 | c task, but the integration was successful as far as t. |

476 | c the problem may have a singularity, or the input |

477 | c may be inappropriate. |

478 | c -5 means there were repeated convergence test failures on |

479 | c one attempted step, before completing the requested |

480 | c task, but the integration was successful as far as t. |

481 | c this may be caused by an inaccurate jacobian matrix, |

482 | c if one is being used. |

483 | c -6 means ewt(i) became zero for some i during the |

484 | c integration. pure relative error control (atol(i)=0.0) |

485 | c was requested on a variable which has now vanished. |

486 | c the integration was successful as far as t. |

487 | c |

488 | c note.. since the normal output value of istate is 2, |

489 | c it does not need to be reset for normal continuation. |

490 | c also, since a negative input value of istate will be |

491 | c regarded as illegal, a negative output value requires the |

492 | c user to change it, and possibly other inputs, before |

493 | c calling the solver again. |

494 | c |

495 | c iopt = an integer flag to specify whether or not any optional |

496 | c inputs are being used on this call. input only. |

497 | c the optional inputs are listed separately below. |

498 | c iopt = 0 means no optional inputs are being used. |

499 | c default values will be used in all cases. |

500 | c iopt = 1 means one or more optional inputs are being used. |

501 | c |

502 | c rwork = a real working array (double precision). |

503 | c the length of rwork must be at least |

504 | c 20 + nyh*(maxord + 1) + 3*neq + lwm where |

505 | c nyh = the initial value of neq, |

506 | c maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a |

507 | c smaller value is given as an optional input), |

508 | c lwm = 0 if miter = 0, |

509 | c lwm = neq**2 + 2 if miter is 1 or 2, |

510 | c lwm = neq + 2 if miter = 3, and |

511 | c lwm = (2*ml+mu+1)*neq + 2 if miter is 4 or 5. |

512 | c (see the mf description for meth and miter.) |

513 | c thus if maxord has its default value and neq is constant, |

514 | c this length is.. |

515 | c 20 + 16*neq for mf = 10, |

516 | c 22 + 16*neq + neq**2 for mf = 11 or 12, |

517 | c 22 + 17*neq for mf = 13, |

518 | c 22 + 17*neq + (2*ml+mu)*neq for mf = 14 or 15, |

519 | c 20 + 9*neq for mf = 20, |

520 | c 22 + 9*neq + neq**2 for mf = 21 or 22, |

521 | c 22 + 10*neq for mf = 23, |

522 | c 22 + 10*neq + (2*ml+mu)*neq for mf = 24 or 25. |

523 | c the first 20 words of rwork are reserved for conditional |

524 | c and optional inputs and optional outputs. |

525 | c |

526 | c the following word in rwork is a conditional input.. |

527 | c rwork(1) = tcrit = critical value of t which the solver |

528 | c is not to overshoot. required if itask is |

529 | c 4 or 5, and ignored otherwise. (see itask.) |

530 | c |

531 | c lrw = the length of the array rwork, as declared by the user. |

532 | c (this will be checked by the solver.) |

533 | c |

534 | c iwork = an integer work array. the length of iwork must be at least |

535 | c 20 if miter = 0 or 3 (mf = 10, 13, 20, 23), or |

536 | c 20 + neq otherwise (mf = 11, 12, 14, 15, 21, 22, 24, 25). |

537 | c the first few words of iwork are used for conditional and |

538 | c optional inputs and optional outputs. |

539 | c |

540 | c the following 2 words in iwork are conditional inputs.. |

541 | c iwork(1) = ml these are the lower and upper |

542 | c iwork(2) = mu half-bandwidths, respectively, of the |

543 | c banded jacobian, excluding the main diagonal. |

544 | c the band is defined by the matrix locations |

545 | c (i,j) with i-ml .le. j .le. i+mu. ml and mu |

546 | c must satisfy 0 .le. ml,mu .le. neq-1. |

547 | c these are required if miter is 4 or 5, and |

548 | c ignored otherwise. ml and mu may in fact be |

549 | c the band parameters for a matrix to which |

550 | c df/dy is only approximately equal. |

551 | c |

552 | c liw = the length of the array iwork, as declared by the user. |

553 | c (this will be checked by the solver.) |

554 | c |

555 | c note.. the work arrays must not be altered between calls to lsode |

556 | c for the same problem, except possibly for the conditional and |

557 | c optional inputs, and except for the last 3*neq words of rwork. |

558 | c the latter space is used for internal scratch space, and so is |

559 | c available for use by the user outside lsode between calls, if |

560 | c desired (but not for use by f or jac). |

561 | c |

562 | c jac = the name of the user-supplied routine (miter = 1 or 4) to |

563 | c compute the jacobian matrix, df/dy, as a function of |

564 | c the scalar t and the vector y. it is to have the form |

565 | c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) |

566 | c dimension y(1), pd(nrowpd,1) |

567 | c where neq, t, y, ml, mu, and nrowpd are input and the array |

568 | c pd is to be loaded with partial derivatives (elements of |

569 | c the jacobian matrix) on output. pd must be given a first |

570 | c dimension of nrowpd. t and y have the same meaning as in |

571 | c subroutine f. (in the dimension statement above, 1 is a |

572 | c dummy dimension.. it can be replaced by any value.) |

573 | c in the full matrix case (miter = 1), ml and mu are |

574 | c ignored, and the jacobian is to be loaded into pd in |

575 | c columnwise manner, with df(i)/dy(j) loaded into pd(i,j). |

576 | c in the band matrix case (miter = 4), the elements |

577 | c within the band are to be loaded into pd in columnwise |

578 | c manner, with diagonal lines of df/dy loaded into the rows |

579 | c of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j). |

580 | c ml and mu are the half-bandwidth parameters (see iwork). |

581 | c the locations in pd in the two triangular areas which |

582 | c correspond to nonexistent matrix elements can be ignored |

583 | c or loaded arbitrarily, as they are overwritten by lsode. |

584 | c jac need not provide df/dy exactly. a crude |

585 | c approximation (possibly with a smaller bandwidth) will do. |

586 | c in either case, pd is preset to zero by the solver, |

587 | c so that only the nonzero elements need be loaded by jac. |

588 | c each call to jac is preceded by a call to f with the same |

589 | c arguments neq, t, and y. thus to gain some efficiency, |

590 | c intermediate quantities shared by both calculations may be |

591 | c saved in a user common block by f and not recomputed by jac, |

592 | c if desired. also, jac may alter the y array, if desired. |

593 | c jac must be declared external in the calling program. |

594 | c subroutine jac may access user-defined quantities in |

595 | c neq(2),... and/or in y(neq(1)+1),... if neq is an array |

596 | c (dimensioned in jac) and/or y has length exceeding neq(1). |

597 | c see the descriptions of neq and y above. |

598 | c |

599 | c mf = the method flag. used only for input. the legal values of |

600 | c mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, and 25. |

601 | c mf has decimal digits meth and miter.. mf = 10*meth + miter. |

602 | c meth indicates the basic linear multistep method.. |

603 | c meth = 1 means the implicit adams method. |

604 | c meth = 2 means the method based on backward |

605 | c differentiation formulas (bdf-s). |

606 | c miter indicates the corrector iteration method.. |

607 | c miter = 0 means functional iteration (no jacobian matrix |

608 | c is involved). |

609 | c miter = 1 means chord iteration with a user-supplied |

610 | c full (neq by neq) jacobian. |

611 | c miter = 2 means chord iteration with an internally |

612 | c generated (difference quotient) full jacobian |

613 | c (using neq extra calls to f per df/dy value). |

614 | c miter = 3 means chord iteration with an internally |

615 | c generated diagonal jacobian approximation. |

616 | c (using 1 extra call to f per df/dy evaluation). |

617 | c miter = 4 means chord iteration with a user-supplied |

618 | c banded jacobian. |

619 | c miter = 5 means chord iteration with an internally |

620 | c generated banded jacobian (using ml+mu+1 extra |

621 | c calls to f per df/dy evaluation). |

622 | c if miter = 1 or 4, the user must supply a subroutine jac |

623 | c (the name is arbitrary) as described above under jac. |

624 | c for other values of miter, a dummy argument can be used. |

625 | c----------------------------------------------------------------------- |

626 | c optional inputs. |

627 | c |

628 | c the following is a list of the optional inputs provided for in the |

629 | c call sequence. (see also part ii.) for each such input variable, |

630 | c this table lists its name as used in this documentation, its |

631 | c location in the call sequence, its meaning, and the default value. |

632 | c the use of any of these inputs requires iopt = 1, and in that |

633 | c case all of these inputs are examined. a value of zero for any |

634 | c of these optional inputs will cause the default value to be used. |

635 | c thus to use a subset of the optional inputs, simply preload |

636 | c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and |

637 | c then set those of interest to nonzero values. |

638 | c |

639 | c name location meaning and default value |

640 | c |

641 | c h0 rwork(5) the step size to be attempted on the first step. |

642 | c the default value is determined by the solver. |

643 | c |

644 | c hmax rwork(6) the maximum absolute step size allowed. |

645 | c the default value is infinite. |

646 | c |

647 | c hmin rwork(7) the minimum absolute step size allowed. |

648 | c the default value is 0. (this lower bound is not |

649 | c enforced on the final step before reaching tcrit |

650 | c when itask = 4 or 5.) |

651 | c |

652 | c maxord iwork(5) the maximum order to be allowed. the default |

653 | c value is 12 if meth = 1, and 5 if meth = 2. |

654 | c if maxord exceeds the default value, it will |

655 | c be reduced to the default value. |

656 | c if maxord is changed during the problem, it may |

657 | c cause the current order to be reduced. |

658 | c |

659 | c mxstep iwork(6) maximum number of (internally defined) steps |

660 | c allowed during one call to the solver. |

661 | c the default value is 500. |

662 | c |

663 | c mxhnil iwork(7) maximum number of messages printed (per problem) |

664 | c warning that t + h = t on a step (h = step size). |

665 | c this must be positive to result in a non-default |

666 | c value. the default value is 10. |

667 | c----------------------------------------------------------------------- |

668 | c optional outputs. |

669 | c |

670 | c as optional additional output from lsode, the variables listed |

671 | c below are quantities related to the performance of lsode |

672 | c which are available to the user. these are communicated by way of |

673 | c the work arrays, but also have internal mnemonic names as shown. |

674 | c except where stated otherwise, all of these outputs are defined |

675 | c on any successful return from lsode, and on any return with |

676 | c istate = -1, -2, -4, -5, or -6. on an illegal input return |

677 | c (istate = -3), they will be unchanged from their existing values |

678 | c (if any), except possibly for tolsf, lenrw, and leniw. |

679 | c on any error return, outputs relevant to the error will be defined, |

680 | c as noted below. |

681 | c |

682 | c name location meaning |

683 | c |

684 | c hu rwork(11) the step size in t last used (successfully). |

685 | c |

686 | c hcur rwork(12) the step size to be attempted on the next step. |

687 | c |

688 | c tcur rwork(13) the current value of the independent variable |

689 | c which the solver has actually reached, i.e. the |

690 | c current internal mesh point in t. on output, tcur |

691 | c will always be at least as far as the argument |

692 | c t, but may be farther (if interpolation was done). |

693 | c |

694 | c tolsf rwork(14) a tolerance scale factor, greater than 1.0, |

695 | c computed when a request for too much accuracy was |

696 | c detected (istate = -3 if detected at the start of |

697 | c the problem, istate = -2 otherwise). if itol is |

698 | c left unaltered but rtol and atol are uniformly |

699 | c scaled up by a factor of tolsf for the next call, |

700 | c then the solver is deemed likely to succeed. |

701 | c (the user may also ignore tolsf and alter the |

702 | c tolerance parameters in any other way appropriate.) |

703 | c |

704 | c nst iwork(11) the number of steps taken for the problem so far. |

705 | c |

706 | c nfe iwork(12) the number of f evaluations for the problem so far. |

707 | c |

708 | c nje iwork(13) the number of jacobian evaluations (and of matrix |

709 | c lu decompositions) for the problem so far. |

710 | c |

711 | c nqu iwork(14) the method order last used (successfully). |

712 | c |

713 | c nqcur iwork(15) the order to be attempted on the next step. |

714 | c |

715 | c imxer iwork(16) the index of the component of largest magnitude in |

716 | c the weighted local error vector ( e(i)/ewt(i) ), |

717 | c on an error return with istate = -4 or -5. |

718 | c |

719 | c lenrw iwork(17) the length of rwork actually required. |

720 | c this is defined on normal returns and on an illegal |

721 | c input return for insufficient storage. |

722 | c |

723 | c leniw iwork(18) the length of iwork actually required. |

724 | c this is defined on normal returns and on an illegal |

725 | c input return for insufficient storage. |

726 | c |

727 | c the following two arrays are segments of the rwork array which |

728 | c may also be of interest to the user as optional outputs. |

729 | c for each array, the table below gives its internal name, |

730 | c its base address in rwork, and its description. |

731 | c |

732 | c name base address description |

733 | c |

734 | c yh 21 the nordsieck history array, of size nyh by |

735 | c (nqcur + 1), where nyh is the initial value |

736 | c of neq. for j = 0,1,...,nqcur, column j+1 |

737 | c of yh contains hcur**j/factorial(j) times |

738 | c the j-th derivative of the interpolating |

739 | c polynomial currently representing the solution, |

740 | c evaluated at t = tcur. |

741 | c |

742 | c acor lenrw-neq+1 array of size neq used for the accumulated |

743 | c corrections on each step, scaled on output |

744 | c to represent the estimated local error in y |

745 | c on the last step. this is the vector e in |

746 | c the description of the error control. it is |

747 | c defined only on a successful return from lsode. |

748 | c |

749 | c----------------------------------------------------------------------- |

750 | c part ii. other routines callable. |

751 | c |

752 | c the following are optional calls which the user may make to |

753 | c gain additional capabilities in conjunction with lsode. |

754 | c (the routines xsetun and xsetf are designed to conform to the |

755 | c slatec error handling package.) |

756 | c |

757 | c form of call function |

758 | c call xsetun(lun) set the logical unit number, lun, for |

759 | c output of messages from lsode, if |

760 | c the default is not desired. |

761 | c the default value of lun is 6. |

762 | c |

763 | c call xsetf(mflag) set a flag to control the printing of |

764 | c messages by lsode. |

765 | c mflag = 0 means do not print. (danger.. |

766 | c this risks losing valuable information.) |

767 | c mflag = 1 means print (the default). |

768 | c |

769 | c either of the above calls may be made at |

770 | c any time and will take effect immediately. |

771 | c |

772 | c call srcom(rsav,isav,job) saves and restores the contents of |

773 | c the internal common blocks used by |

774 | c lsode (see part iii below). |

775 | c rsav must be a real array of length 218 |

776 | c or more, and isav must be an integer |

777 | c array of length 41 or more. |

778 | c job=1 means save common into rsav/isav. |

779 | c job=2 means restore common from rsav/isav. |

780 | c srcom is useful if one is |

781 | c interrupting a run and restarting |

782 | c later, or alternating between two or |

783 | c more problems solved with lsode. |

784 | c |

785 | c call intdy(,,,,,) provide derivatives of y, of various |

786 | c (see below) orders, at a specified point t, if |

787 | c desired. it may be called only after |

788 | c a successful return from lsode. |

789 | c |

790 | c the detailed instructions for using intdy are as follows. |

791 | c the form of the call is.. |

792 | c |

793 | c call intdy (t, k, rwork(21), nyh, dky, iflag) |

794 | c |

795 | c the input parameters are.. |

796 | c |

797 | c t = value of independent variable where answers are desired |

798 | c (normally the same as the t last returned by lsode). |

799 | c for valid results, t must lie between tcur - hu and tcur. |

800 | c (see optional outputs for tcur and hu.) |

801 | c k = integer order of the derivative desired. k must satisfy |

802 | c 0 .le. k .le. nqcur, where nqcur is the current order |

803 | c (see optional outputs). the capability corresponding |

804 | c to k = 0, i.e. computing y(t), is already provided |

805 | c by lsode directly. since nqcur .ge. 1, the first |

806 | c derivative dy/dt is always available with intdy. |

807 | c rwork(21) = the base address of the history array yh. |

808 | c nyh = column length of yh, equal to the initial value of neq. |

809 | c |

810 | c the output parameters are.. |

811 | c |

812 | c dky = a real array of length neq containing the computed value |

813 | c of the k-th derivative of y(t). |

814 | c iflag = integer flag, returned as 0 if k and t were legal, |

815 | c -1 if k was illegal, and -2 if t was illegal. |

816 | c on an error return, a message is also written. |

817 | c----------------------------------------------------------------------- |

818 | c part iii. common blocks. |

819 | c |

820 | c if lsode is to be used in an overlay situation, the user |

821 | c must declare, in the primary overlay, the variables in.. |

822 | c (1) the call sequence to lsode, |

823 | c (2) the two internal common blocks |

824 | c /ls0001/ of length 257 (218 double precision words |

825 | c followed by 39 integer words), |

826 | c /eh0001/ of length 2 (integer words). |

827 | c |

828 | c if lsode is used on a system in which the contents of internal |

829 | c common blocks are not preserved between calls, the user should |

830 | c declare the above two common blocks in his main program to insure |

831 | c that their contents are preserved. |

832 | c |

833 | c if the solution of a given problem by lsode is to be interrupted |

834 | c and then later continued, such as when restarting an interrupted run |

835 | c or alternating between two or more problems, the user should save, |

836 | c following the return from the last lsode call prior to the |

837 | c interruption, the contents of the call sequence variables and the |

838 | c internal common blocks, and later restore these values before the |

839 | c next lsode call for that problem. to save and restore the common |

840 | c blocks, use subroutine srcom (see part ii above). |

841 | c |

842 | c----------------------------------------------------------------------- |

843 | c part iv. optionally replaceable solver routines. |

844 | c |

845 | c below are descriptions of two routines in the lsode package which |

846 | c relate to the measurement of errors. either routine can be |

847 | c replaced by a user-supplied version, if desired. however, since such |

848 | c a replacement may have a major impact on performance, it should be |

849 | c done only when absolutely necessary, and only with great caution. |

850 | c (note.. the means by which the package version of a routine is |

851 | c superseded by the user-s version may be system-dependent.) |

852 | c |

853 | c (a) ewset. |

854 | c the following subroutine is called just before each internal |

855 | c integration step, and sets the array of error weights, ewt, as |

856 | c described under itol/rtol/atol above.. |

857 | c subroutine ewset (neq, itol, rtol, atol, ycur, ewt) |

858 | c where neq, itol, rtol, and atol are as in the lsode call sequence, |

859 | c ycur contains the current dependent variable vector, and |

860 | c ewt is the array of weights set by ewset. |

861 | c |

862 | c if the user supplies this subroutine, it must return in ewt(i) |

863 | c (i = 1,...,neq) a positive quantity suitable for comparing errors |

864 | c in y(i) to. the ewt array returned by ewset is passed to the |

865 | c vnorm routine (see below), and also used by lsode in the computation |

866 | c of the optional output imxer, the diagonal jacobian approximation, |

867 | c and the increments for difference quotient jacobians. |

868 | c |

869 | c in the user-supplied version of ewset, it may be desirable to use |

870 | c the current values of derivatives of y. derivatives up to order nq |

871 | c are available from the history array yh, described above under |

872 | c optional outputs. in ewset, yh is identical to the ycur array, |

873 | c extended to nq + 1 columns with a column length of nyh and scale |

874 | c factors of h**j/factorial(j). on the first call for the problem, |

875 | c given by nst = 0, nq is 1 and h is temporarily set to 1.0. |

876 | c the quantities nq, nyh, h, and nst can be obtained by including |

877 | c in ewset the statements.. |

878 | c double precision h, rls |

879 | c common /ls0001/ rls(218),ils(39) |

880 | c nq = ils(35) |

881 | c nyh = ils(14) |

882 | c nst = ils(36) |

883 | c h = rls(212) |

884 | c thus, for example, the current value of dy/dt can be obtained as |

885 | c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is |

886 | c unnecessary when nst = 0). |

887 | c |

888 | c (b) vnorm. |

889 | c the following is a real function routine which computes the weighted |

890 | c root-mean-square norm of a vector v.. |

891 | c d = vnorm (n, v, w) |

892 | c where.. |

893 | c n = the length of the vector, |

894 | c v = real array of length n containing the vector, |

895 | c w = real array of length n containing weights, |

896 | c d = sqrt( (1/n) * sum(v(i)*w(i))**2 ). |

897 | c vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where |

898 | c ewt is as set by subroutine ewset. |

899 | c |

900 | c if the user supplies this function, it should return a non-negative |

901 | c value of vnorm suitable for use in the error control in lsode. |

902 | c none of the arguments should be altered by vnorm. |

903 | c for example, a user-supplied vnorm routine might.. |

904 | c -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or |

905 | c -ignore some components of v in the norm, with the effect of |

906 | c suppressing the error control on those components of y. |

907 | c----------------------------------------------------------------------- |

908 | c----------------------------------------------------------------------- |

909 | c other routines in the lsode package. |

910 | c |

911 | c in addition to subroutine lsode, the lsode package includes the |

912 | c following subroutines and function routines.. |

913 | c intdy computes an interpolated value of the y vector at t = tout. |

914 | c stode is the core integrator, which does one step of the |

915 | c integration and the associated error control. |

916 | c cfode sets all method coefficients and test constants. |

917 | c prepj computes and preprocesses the jacobian matrix j = df/dy |

918 | c and the newton iteration matrix p = i - h*l0*j. |

919 | c solsy manages solution of linear system in chord iteration. |

920 | c ewset sets the error weight vector ewt before each step. |

921 | c vnorm computes the weighted r.m.s. norm of a vector. |

922 | c srcom is a user-callable routine to save and restore |

923 | c the contents of the internal common blocks. |

924 | c dgefa and dgesl are routines from linpack for solving full |

925 | c systems of linear algebraic equations. |

926 | c dgbfa and dgbsl are routines from linpack for solving banded |

927 | c linear systems. |

928 | c daxpy, dscal, idamax, and ddot are basic linear algebra modules |

929 | c (blas) used by the above linpack routines. |

930 | c d1mach computes the unit roundoff in a machine-independent manner. |

931 | c XERRWV, xsetun, and xsetf handle the printing of all error |

932 | c messages and warnings. XERRWV is machine-dependent. |

933 | C ASCEND C XERRWV has been replaced by XASCWV from interface/Lsode.c. |

934 | c note.. vnorm, idamax, ddot, and d1mach are function routines. |

935 | c all the others are subroutines. |

936 | c |

937 | c the intrinsic and external routines used by lsode are.. |

938 | c dabs, dmax1, dmin1, dfloat, max0, min0, mod, dsign, dsqrt, and write. |

939 | c |

940 | c a block data subprogram is also included with the package, |

941 | c for loading some of the variables in internal common. |

942 | c |

943 | c----------------------------------------------------------------------- |

944 | c the following card is for optimized compilation on llnl compilers. |

945 | clll. optimize |

946 | c----------------------------------------------------------------------- |

947 | cascend changes |

948 | Ckaa external aftime |

949 | external prepj, solsy |

950 | integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, |

951 | 1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns |

952 | integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

953 | 1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

954 | integer i, i1, i2, iflag, imxer, kgo, lf0, |

955 | 1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0 |

956 | double precision rowns, |

957 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |

958 | double precision atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli, |

959 | 1 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0, |

960 | 2 d1mach, vnorm |

961 | dimension mord(2) |

962 | logical ihit |

963 | c----------------------------------------------------------------------- |

964 | c the following internal common block contains |

965 | c (a) variables which are local to any subroutine but whose values must |

966 | c be preserved between calls to the routine (own variables), and |

967 | c (b) variables which are communicated between subroutines. |

968 | c the structure of the block is as follows.. all real variables are |

969 | c listed first, followed by all integers. within each type, the |

970 | c variables are grouped with those local to subroutine lsode first, |

971 | c then those local to subroutine stode, and finally those used |

972 | c for communication. the block is declared in subroutines |

973 | c lsode, intdy, stode, prepj, and solsy. groups of variables are |

974 | c replaced by dummy arrays in the common declarations in routines |

975 | c where those variables are not used. |

976 | c----------------------------------------------------------------------- |

977 | common /ls0001/ rowns(209), |

978 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |

979 | 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, |

980 | 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), |

981 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

982 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

983 | c |

984 | data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/ |

985 | |

986 | C ascend debugging options |

987 | C type *, 'iopt',iopt |

988 | C type *, 'h0',rwork(5) |

989 | C type *, 'hmax',rwork(6) |

990 | C type *, 'hmin',rwork(7) |

991 | C type *, 'maxs',iwork(6) |

992 | c----------------------------------------------------------------------- |

993 | c block a. |

994 | c this code block is executed on every call. |

995 | c it tests istate and itask for legality and branches appropriately. |

996 | c if istate .gt. 1 but the flag init shows that initialization has |

997 | c not yet been done, an error return occurs. |

998 | c if istate = 1 and tout = t, jump to block g and return immediately. |

999 | c----------------------------------------------------------------------- |

1000 | if (istate .lt. 1 .or. istate .gt. 3) go to 601 |

1001 | if (itask .lt. 1 .or. itask .gt. 5) go to 602 |

1002 | if (istate .eq. 1) go to 10 |

1003 | if (init .eq. 0) go to 603 |

1004 | if (istate .eq. 2) go to 200 |

1005 | go to 20 |

1006 | 10 init = 0 |

1007 | if (tout .eq. t) go to 430 |

1008 | 20 ntrep = 0 |

1009 | c----------------------------------------------------------------------- |

1010 | c block b. |

1011 | c the next code block is executed for the initial call (istate = 1), |

1012 | c or for a continuation call with parameter changes (istate = 3). |

1013 | c it contains checking of all inputs and various initializations. |

1014 | c |

1015 | c first check legality of the non-optional inputs neq, itol, iopt, |

1016 | c mf, ml, and mu. |

1017 | c----------------------------------------------------------------------- |

1018 | if (neq(1) .le. 0) go to 604 |

1019 | if (istate .eq. 1) go to 25 |

1020 | if (neq(1) .gt. n) go to 605 |

1021 | 25 n = neq(1) |

1022 | if (itol .lt. 1 .or. itol .gt. 4) go to 606 |

1023 | if (iopt .lt. 0 .or. iopt .gt. 1) go to 607 |

1024 | meth = mf/10 |

1025 | miter = mf - 10*meth |

1026 | if (meth .lt. 1 .or. meth .gt. 2) go to 608 |

1027 | if (miter .lt. 0 .or. miter .gt. 5) go to 608 |

1028 | if (miter .le. 3) go to 30 |

1029 | ml = iwork(1) |

1030 | mu = iwork(2) |

1031 | if (ml .lt. 0 .or. ml .ge. n) go to 609 |

1032 | if (mu .lt. 0 .or. mu .ge. n) go to 610 |

1033 | 30 continue |

1034 | c next process and check the optional inputs. -------------------------- |

1035 | if (iopt .eq. 1) go to 40 |

1036 | maxord = mord(meth) |

1037 | mxstep = mxstp0 |

1038 | mxhnil = mxhnl0 |

1039 | if (istate .eq. 1) h0 = 0.0d0 |

1040 | hmxi = 0.0d0 |

1041 | hmin = 0.0d0 |

1042 | go to 60 |

1043 | 40 maxord = iwork(5) |

1044 | if (maxord .lt. 0) go to 611 |

1045 | if (maxord .eq. 0) maxord = 100 |

1046 | maxord = min0(maxord,mord(meth)) |

1047 | mxstep = iwork(6) |

1048 | if (mxstep .lt. 0) go to 612 |

1049 | if (mxstep .eq. 0) mxstep = mxstp0 |

1050 | mxhnil = iwork(7) |

1051 | if (mxhnil .lt. 0) go to 613 |

1052 | if (mxhnil .eq. 0) mxhnil = mxhnl0 |

1053 | if (istate .ne. 1) go to 50 |

1054 | h0 = rwork(5) |

1055 | if ((tout - t)*h0 .lt. 0.0d0) go to 614 |

1056 | 50 hmax = rwork(6) |

1057 | if (hmax .lt. 0.0d0) go to 615 |

1058 | hmxi = 0.0d0 |

1059 | if (hmax .gt. 0.0d0) hmxi = 1.0d0/hmax |

1060 | hmin = rwork(7) |

1061 | if (hmin .lt. 0.0d0) go to 616 |

1062 | c----------------------------------------------------------------------- |

1063 | c set work array pointers and check lengths lrw and liw. |

1064 | c pointers to segments of rwork and iwork are named by prefixing l to |

1065 | c the name of the segment. e.g., the segment yh starts at rwork(lyh). |

1066 | c segments of rwork (in order) are denoted yh, wm, ewt, savf, acor. |

1067 | c----------------------------------------------------------------------- |

1068 | 60 lyh = 21 |

1069 | if (istate .eq. 1) nyh = n |

1070 | lwm = lyh + (maxord + 1)*nyh |

1071 | if (miter .eq. 0) lenwm = 0 |

1072 | if (miter .eq. 1 .or. miter .eq. 2) lenwm = n*n + 2 |

1073 | if (miter .eq. 3) lenwm = n + 2 |

1074 | if (miter .ge. 4) lenwm = (2*ml + mu + 1)*n + 2 |

1075 | lewt = lwm + lenwm |

1076 | lsavf = lewt + n |

1077 | lacor = lsavf + n |

1078 | lenrw = lacor + n - 1 |

1079 | iwork(17) = lenrw |

1080 | liwm = 1 |

1081 | leniw = 20 + n |

1082 | if (miter .eq. 0 .or. miter .eq. 3) leniw = 20 |

1083 | iwork(18) = leniw |

1084 | if (lenrw .gt. lrw) go to 617 |

1085 | if (leniw .gt. liw) go to 618 |

1086 | c check rtol and atol for legality. ------------------------------------ |

1087 | rtoli = rtol(1) |

1088 | atoli = atol(1) |

1089 | do 70 i = 1,n |

1090 | if (itol .ge. 3) rtoli = rtol(i) |

1091 | if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) |

1092 | if (rtoli .lt. 0.0d0) go to 619 |

1093 | if (atoli .lt. 0.0d0) go to 620 |

1094 | 70 continue |

1095 | if (istate .eq. 1) go to 100 |

1096 | c if istate = 3, set flag to signal parameter changes to stode. -------- |

1097 | jstart = -1 |

1098 | if (nq .le. maxord) go to 90 |

1099 | c maxord was reduced below nq. copy yh(*,maxord+2) into savf. --------- |

1100 | do 80 i = 1,n |

1101 | 80 rwork(i+lsavf-1) = rwork(i+lwm-1) |

1102 | c reload wm(1) = rwork(lwm), since lwm may have changed. --------------- |

1103 | 90 if (miter .gt. 0) rwork(lwm) = dsqrt(uround) |

1104 | if (n .eq. nyh) go to 200 |

1105 | c neq was reduced. zero part of yh to avoid undefined references. ----- |

1106 | i1 = lyh + l*nyh |

1107 | i2 = lyh + (maxord + 1)*nyh - 1 |

1108 | if (i1 .gt. i2) go to 200 |

1109 | do 95 i = i1,i2 |

1110 | 95 rwork(i) = 0.0d0 |

1111 | go to 200 |

1112 | c----------------------------------------------------------------------- |

1113 | c block c. |

1114 | c the next block is for the initial call only (istate = 1). |

1115 | c it contains all remaining initializations, the initial call to f, |

1116 | c and the calculation of the initial step size. |

1117 | c the error weights in ewt are inverted after being loaded. |

1118 | c----------------------------------------------------------------------- |

1119 | 100 uround = d1mach(4) |

1120 | tn = t |

1121 | if (itask .ne. 4 .and. itask .ne. 5) go to 110 |

1122 | tcrit = rwork(1) |

1123 | if ((tcrit - tout)*(tout - t) .lt. 0.0d0) go to 625 |

1124 | if (h0 .ne. 0.0d0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d0) |

1125 | 1 h0 = tcrit - t |

1126 | 110 jstart = 0 |

1127 | if (miter .gt. 0) rwork(lwm) = dsqrt(uround) |

1128 | nhnil = 0 |

1129 | nst = 0 |

1130 | nje = 0 |

1131 | nslast = 0 |

1132 | hu = 0.0d0 |

1133 | nqu = 0 |

1134 | ccmax = 0.3d0 |

1135 | maxcor = 3 |

1136 | msbp = 20 |

1137 | mxncf = 10 |

1138 | c initial call to f. (lf0 points to yh(*,2).) ------------------------- |

1139 | lf0 = lyh + nyh |

1140 | call f (neq, t, y, rwork(lf0)) |

1141 | nfe = 1 |

1142 | c load the initial value vector in yh. --------------------------------- |

1143 | do 115 i = 1,n |

1144 | 115 rwork(i+lyh-1) = y(i) |

1145 | c load and invert the ewt array. (h is temporarily set to 1.0.) ------- |

1146 | nq = 1 |

1147 | h = 1.0d0 |

1148 | call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) |

1149 | do 120 i = 1,n |

1150 | if (rwork(i+lewt-1) .le. 0.0d0) go to 621 |

1151 | 120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1) |

1152 | c----------------------------------------------------------------------- |

1153 | c the coding below computes the step size, h0, to be attempted on the |

1154 | c first step, unless the user has supplied a value for this. |

1155 | c first check that tout - t differs significantly from zero. |

1156 | c a scalar tolerance quantity tol is computed, as max(rtol(i)) |

1157 | c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted |

1158 | c so as to be between 100*uround and 1.0e-3. |

1159 | c then the computed value h0 is given by.. |

1160 | c neq |

1161 | c h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2 ) |

1162 | c 1 |

1163 | c where w0 = max ( abs(t), abs(tout) ), |

1164 | c f(i) = i-th component of initial value of f, |

1165 | c ywt(i) = ewt(i)/tol (a weight for y(i)). |

1166 | c the sign of h0 is inferred from the initial values of tout and t. |

1167 | c----------------------------------------------------------------------- |

1168 | if (h0 .ne. 0.0d0) go to 180 |

1169 | tdist = dabs(tout - t) |

1170 | w0 = dmax1(dabs(t),dabs(tout)) |

1171 | if (tdist .lt. 2.0d0*uround*w0) go to 622 |

1172 | tol = rtol(1) |

1173 | if (itol .le. 2) go to 140 |

1174 | do 130 i = 1,n |

1175 | 130 tol = dmax1(tol,rtol(i)) |

1176 | 140 if (tol .gt. 0.0d0) go to 160 |

1177 | atoli = atol(1) |

1178 | do 150 i = 1,n |

1179 | if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) |

1180 | ayi = dabs(y(i)) |

1181 | if (ayi .ne. 0.0d0) tol = dmax1(tol,atoli/ayi) |

1182 | 150 continue |

1183 | 160 tol = dmax1(tol,100.0d0*uround) |

1184 | tol = dmin1(tol,0.001d0) |

1185 | sum = vnorm (n, rwork(lf0), rwork(lewt)) |

1186 | sum = 1.0d0/(tol*w0*w0) + tol*sum**2 |

1187 | h0 = 1.0d0/dsqrt(sum) |

1188 | h0 = dmin1(h0,tdist) |

1189 | h0 = dsign(h0,tout-t) |

1190 | c adjust h0 if necessary to meet hmax bound. --------------------------- |

1191 | 180 rh = dabs(h0)*hmxi |

1192 | if (rh .gt. 1.0d0) h0 = h0/rh |

1193 | c load h with h0 and scale yh(*,2) by h0. ------------------------------ |

1194 | h = h0 |

1195 | do 190 i = 1,n |

1196 | 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1) |

1197 | go to 270 |

1198 | c----------------------------------------------------------------------- |

1199 | c block d. |

1200 | c the next code block is for continuation calls only (istate = 2 or 3) |

1201 | c and is to check stop conditions before taking a step. |

1202 | c----------------------------------------------------------------------- |

1203 | 200 nslast = nst |

1204 | go to (210, 250, 220, 230, 240), itask |

1205 | 210 if ((tn - tout)*h .lt. 0.0d0) go to 250 |

1206 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |

1207 | if (iflag .ne. 0) go to 627 |

1208 | t = tout |

1209 | go to 420 |

1210 | 220 tp = tn - hu*(1.0d0 + 100.0d0*uround) |

1211 | if ((tp - tout)*h .gt. 0.0d0) go to 623 |

1212 | if ((tn - tout)*h .lt. 0.0d0) go to 250 |

1213 | go to 400 |

1214 | 230 tcrit = rwork(1) |

1215 | if ((tn - tcrit)*h .gt. 0.0d0) go to 624 |

1216 | if ((tcrit - tout)*h .lt. 0.0d0) go to 625 |

1217 | if ((tn - tout)*h .lt. 0.0d0) go to 245 |

1218 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |

1219 | if (iflag .ne. 0) go to 627 |

1220 | t = tout |

1221 | go to 420 |

1222 | 240 tcrit = rwork(1) |

1223 | if ((tn - tcrit)*h .gt. 0.0d0) go to 624 |

1224 | 245 hmx = dabs(tn) + dabs(h) |

1225 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |

1226 | if (ihit) go to 400 |

1227 | tnext = tn + h*(1.0d0 + 4.0d0*uround) |

1228 | if ((tnext - tcrit)*h .le. 0.0d0) go to 250 |

1229 | h = (tcrit - tn)*(1.0d0 - 4.0d0*uround) |

1230 | if (istate .eq. 2) jstart = -2 |

1231 | c----------------------------------------------------------------------- |

1232 | c block e. |

1233 | c the next block is normally executed for all calls and contains |

1234 | c the call to the one-step core integrator stode. |

1235 | c |

1236 | c this is a looping point for the integration steps. |

1237 | c |

1238 | c first check for too many steps being taken, update ewt (if not at |

1239 | c start of problem), check for too much accuracy being requested, and |

1240 | c check for h below the roundoff level in t. |

1241 | c----------------------------------------------------------------------- |

1242 | 250 continue |

1243 | if ((nst-nslast) .ge. mxstep) go to 500 |

1244 | call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) |

1245 | do 260 i = 1,n |

1246 | if (rwork(i+lewt-1) .le. 0.0d0) go to 510 |

1247 | 260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1) |

1248 | 270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt)) |

1249 | if (tolsf .le. 1.0d0) go to 280 |

1250 | tolsf = tolsf*2.0d0 |

1251 | if (nst .eq. 0) go to 626 |

1252 | go to 520 |

1253 | 280 if ((tn + h) .ne. tn) go to 290 |

1254 | nhnil = nhnil + 1 |

1255 | if (nhnil .gt. mxhnil) go to 290 |

1256 | call xascwv(50hlsode-- warning..internal t (=r1) and h (=r2) are, |

1257 | 1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1258 | call xascwv( |

1259 | 1 60h such that in the machine, t + h = t on the next step , |

1260 | 1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1261 | call xascwv(50h (h = step size). solver will continue anyway, |

1262 | 1 50, 101, 0, 0, 0, 0, 2, tn, h) |

1263 | if (nhnil .lt. mxhnil) go to 290 |

1264 | call xascwv(50hlsode-- above warning has been issued i1 times. , |

1265 | 1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1266 | call xascwv(50h it will not be issued again for this problem, |

1267 | 1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0) |

1268 | 290 continue |

1269 | c----------------------------------------------------------------------- |

1270 | c call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,jac,prepj,solsy) |

1271 | c----------------------------------------------------------------------- |

1272 | c boyd |

1273 | call stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt), |

1274 | 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm), |

1275 | 2 f, jac, prepj, solsy) |

1276 | kgo = 1 - kflag |

1277 | go to (300, 530, 540), kgo |

1278 | c----------------------------------------------------------------------- |

1279 | c block f. |

1280 | c the following block handles the case of a successful return from the |

1281 | c core integrator (kflag = 0). test for stop conditions. |

1282 | c----------------------------------------------------------------------- |

1283 | 300 init = 1 |

1284 | go to (310, 400, 330, 340, 350), itask |

1285 | c itask = 1. if tout has been reached, interpolate. ------------------- |

1286 | 310 if ((tn - tout)*h .lt. 0.0d0) go to 250 |

1287 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |

1288 | t = tout |

1289 | go to 420 |

1290 | c itask = 3. jump to exit if tout was reached. ------------------------ |

1291 | 330 if ((tn - tout)*h .ge. 0.0d0) go to 400 |

1292 | go to 250 |

1293 | c itask = 4. see if tout or tcrit was reached. adjust h if necessary. |

1294 | 340 if ((tn - tout)*h .lt. 0.0d0) go to 345 |

1295 | call intdy (tout, 0, rwork(lyh), nyh, y, iflag) |

1296 | t = tout |

1297 | go to 420 |

1298 | 345 hmx = dabs(tn) + dabs(h) |

1299 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |

1300 | if (ihit) go to 400 |

1301 | tnext = tn + h*(1.0d0 + 4.0d0*uround) |

1302 | if ((tnext - tcrit)*h .le. 0.0d0) go to 250 |

1303 | h = (tcrit - tn)*(1.0d0 - 4.0d0*uround) |

1304 | jstart = -2 |

1305 | go to 250 |

1306 | c itask = 5. see if tcrit was reached and jump to exit. --------------- |

1307 | 350 hmx = dabs(tn) + dabs(h) |

1308 | ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx |

1309 | c----------------------------------------------------------------------- |

1310 | c block g. |

1311 | c the following block handles all successful returns from lsode. |

1312 | c if itask .ne. 1, y is loaded from yh and t is set accordingly. |

1313 | c istate is set to 2, the illegal input counter is zeroed, and the |

1314 | c optional outputs are loaded into the work arrays before returning. |

1315 | c if istate = 1 and tout = t, there is a return with no action taken, |

1316 | c except that if this has happened repeatedly, the run is terminated. |

1317 | c----------------------------------------------------------------------- |

1318 | 400 do 410 i = 1,n |

1319 | 410 y(i) = rwork(i+lyh-1) |

1320 | t = tn |

1321 | if (itask .ne. 4 .and. itask .ne. 5) go to 420 |

1322 | if (ihit) t = tcrit |

1323 | 420 istate = 2 |

1324 | illin = 0 |

1325 | rwork(11) = hu |

1326 | rwork(12) = h |

1327 | rwork(13) = tn |

1328 | iwork(11) = nst |

1329 | iwork(12) = nfe |

1330 | iwork(13) = nje |

1331 | iwork(14) = nqu |

1332 | iwork(15) = nq |

1333 | return |

1334 | c |

1335 | 430 ntrep = ntrep + 1 |

1336 | if (ntrep .lt. 5) return |

1337 | call xascwv( |

1338 | 1 60hlsode-- repeated calls with istate = 1 and tout = t (=r1) , |

1339 | 1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0) |

1340 | go to 800 |

1341 | c----------------------------------------------------------------------- |

1342 | c block h. |

1343 | c the following block handles all unsuccessful returns other than |

1344 | c those for illegal input. first the error message routine is called. |

1345 | c if there was an error test or convergence test failure, imxer is set. |

1346 | c then y is loaded from yh, t is set to tn, and the illegal input |

1347 | c counter illin is set to 0. the optional outputs are loaded into |

1348 | c the work arrays before returning. |

1349 | c----------------------------------------------------------------------- |

1350 | c the maximum number of steps was taken before reaching tout. ---------- |

1351 | 500 call xascwv(50hlsode-- at current t (=r1), mxstep (=i1) steps , |

1352 | 1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1353 | call xascwv(50h taken on this call before reaching tout , |

1354 | 1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0) |

1355 | istate = -1 |

1356 | go to 580 |

1357 | c ewt(i) .le. 0.0 for some i (not at start of problem). ---------------- |

1358 | 510 ewti = rwork(lewt+i-1) |

1359 | call xascwv(50hlsode-- at t (=r1), ewt(i1) has become r2 .le. 0., |

1360 | 1 50, 202, 0, 1, i, 0, 2, tn, ewti) |

1361 | istate = -6 |

1362 | go to 580 |

1363 | c too much accuracy requested for machine precision. ------------------- |

1364 | 520 call xascwv(50hlsode-- at t (=r1), too much accuracy requested , |

1365 | 1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1366 | call xascwv(50h for precision of machine.. see tolsf (=r2) , |

1367 | 1 50, 203, 0, 0, 0, 0, 2, tn, tolsf) |

1368 | rwork(14) = tolsf |

1369 | istate = -2 |

1370 | go to 580 |

1371 | c kflag = -1. error test failed repeatedly or with abs(h) = hmin. ----- |

1372 | 530 call xascwv(50hlsode-- at t(=r1) and step size h(=r2), the error, |

1373 | 1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1374 | call xascwv(50h test failed repeatedly or with abs(h) = hmin, |

1375 | 1 50, 204, 0, 0, 0, 0, 2, tn, h) |

1376 | istate = -4 |

1377 | go to 560 |

1378 | c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ---- |

1379 | 540 call xascwv(50hlsode-- at t (=r1) and step size h (=r2), the , |

1380 | 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1381 | call xascwv(50h corrector convergence failed repeatedly , |

1382 | 1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1383 | call xascwv(30h or with abs(h) = hmin , |

1384 | 1 30, 205, 0, 0, 0, 0, 2, tn, h) |

1385 | istate = -5 |

1386 | c compute imxer if relevant. ------------------------------------------- |

1387 | 560 big = 0.0d0 |

1388 | imxer = 1 |

1389 | do 570 i = 1,n |

1390 | size = dabs(rwork(i+lacor-1)*rwork(i+lewt-1)) |

1391 | if (big .ge. size) go to 570 |

1392 | big = size |

1393 | imxer = i |

1394 | 570 continue |

1395 | iwork(16) = imxer |

1396 | c set y vector, t, illin, and optional outputs. ------------------------ |

1397 | 580 do 590 i = 1,n |

1398 | 590 y(i) = rwork(i+lyh-1) |

1399 | t = tn |

1400 | illin = 0 |

1401 | rwork(11) = hu |

1402 | rwork(12) = h |

1403 | rwork(13) = tn |

1404 | iwork(11) = nst |

1405 | iwork(12) = nfe |

1406 | iwork(13) = nje |

1407 | iwork(14) = nqu |

1408 | iwork(15) = nq |

1409 | return |

1410 | c----------------------------------------------------------------------- |

1411 | c block i. |

1412 | c the following block handles all error returns due to illegal input |

1413 | c (istate = -3), as detected before calling the core integrator. |

1414 | c first the error message routine is called. then if there have been |

1415 | c 5 consecutive such returns just before this call to the solver, |

1416 | c the run is halted. |

1417 | c----------------------------------------------------------------------- |

1418 | 601 call xascwv(30hlsode-- istate (=i1) illegal , |

1419 | 1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0) |

1420 | go to 700 |

1421 | 602 call xascwv(30hlsode-- itask (=i1) illegal , |

1422 | 1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0) |

1423 | go to 700 |

1424 | 603 call xascwv(50hlsode-- istate .gt. 1 but lsode not initialized , |

1425 | 1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1426 | go to 700 |

1427 | 604 call xascwv(30hlsode-- neq (=i1) .lt. 1 , |

1428 | 1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0) |

1429 | go to 700 |

1430 | 605 call xascwv(50hlsode-- istate = 3 and neq increased (i1 to i2) , |

1431 | 1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0) |

1432 | go to 700 |

1433 | 606 call xascwv(30hlsode-- itol (=i1) illegal , |

1434 | 1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0) |

1435 | go to 700 |

1436 | 607 call xascwv(30hlsode-- iopt (=i1) illegal , |

1437 | 1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0) |

1438 | go to 700 |

1439 | 608 call xascwv(30hlsode-- mf (=i1) illegal , |

1440 | 1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0) |

1441 | go to 700 |

1442 | 609 call xascwv(50hlsode-- ml (=i1) illegal.. .lt.0 or .ge.neq (=i2), |

1443 | 1 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0) |

1444 | go to 700 |

1445 | 610 call xascwv(50hlsode-- mu (=i1) illegal.. .lt.0 or .ge.neq (=i2), |

1446 | 1 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0) |

1447 | go to 700 |

1448 | 611 call xascwv(30hlsode-- maxord (=i1) .lt. 0 , |

1449 | 1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0) |

1450 | go to 700 |

1451 | 612 call xascwv(30hlsode-- mxstep (=i1) .lt. 0 , |

1452 | 1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0) |

1453 | go to 700 |

1454 | 613 call xascwv(30hlsode-- mxhnil (=i1) .lt. 0 , |

1455 | 1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0) |

1456 | go to 700 |

1457 | 614 call xascwv(40hlsode-- tout (=r1) behind t (=r2) , |

1458 | 1 40, 14, 0, 0, 0, 0, 2, tout, t) |

1459 | call xascwv(50h integration direction is given by h0 (=r1) , |

1460 | 1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0) |

1461 | go to 700 |

1462 | 615 call xascwv(30hlsode-- hmax (=r1) .lt. 0.0 , |

1463 | 1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0) |

1464 | go to 700 |

1465 | 616 call xascwv(30hlsode-- hmin (=r1) .lt. 0.0 , |

1466 | 1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0) |

1467 | go to 700 |

1468 | 617 call xascwv( |

1469 | 1 60hlsode-- rwork length needed, lenrw (=i1), exceeds lrw (=i2), |

1470 | 1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0) |

1471 | go to 700 |

1472 | 618 call xascwv( |

1473 | 1 60hlsode-- iwork length needed, leniw (=i1), exceeds liw (=i2), |

1474 | 1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0) |

1475 | go to 700 |

1476 | 619 call xascwv(40hlsode-- rtol(i1) is r1 .lt. 0.0 , |

1477 | 1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0) |

1478 | go to 700 |

1479 | 620 call xascwv(40hlsode-- atol(i1) is r1 .lt. 0.0 , |

1480 | 1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0) |

1481 | go to 700 |

1482 | 621 ewti = rwork(lewt+i-1) |

1483 | call xascwv(40hlsode-- ewt(i1) is r1 .le. 0.0 , |

1484 | 1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0) |

1485 | go to 700 |

1486 | 622 call xascwv( |

1487 | 1 60hlsode-- tout (=r1) too close to t(=r2) to start integration, |

1488 | 1 60, 22, 0, 0, 0, 0, 2, tout, t) |

1489 | go to 700 |

1490 | 623 call xascwv( |

1491 | 1 60hlsode-- itask = i1 and tout (=r1) behind tcur - hu (= r2) , |

1492 | 1 60, 23, 0, 1, itask, 0, 2, tout, tp) |

1493 | go to 700 |

1494 | 624 call xascwv( |

1495 | 1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) , |

1496 | 1 60, 24, 0, 0, 0, 0, 2, tcrit, tn) |

1497 | go to 700 |

1498 | 625 call xascwv( |

1499 | 1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) , |

1500 | 1 60, 25, 0, 0, 0, 0, 2, tcrit, tout) |

1501 | go to 700 |

1502 | 626 call xascwv(50hlsode-- at start of problem, too much accuracy , |

1503 | 1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1504 | call xascwv( |

1505 | 1 60h requested for precision of machine.. see tolsf (=r1) , |

1506 | 1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0) |

1507 | rwork(14) = tolsf |

1508 | go to 700 |

1509 | 627 call xascwv(50hlsode-- trouble from intdy. itask = i1, tout = r1, |

1510 | 1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0) |

1511 | c |

1512 | 700 if (illin .eq. 5) go to 710 |

1513 | illin = illin + 1 |

1514 | istate = -3 |

1515 | return |

1516 | 710 call xascwv(50hlsode-- repeated occurrences of illegal input , |

1517 | 1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1518 | c |

1519 | 800 call xascwv(50hlsode-- run aborted.. apparent infinite loop , |

1520 | 1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0) |

1521 | return |

1522 | c----------------------- end of subroutine lsode ----------------------- |

1523 | end |

1524 | subroutine solsy (wm, iwm, x, tem) |

1525 | clll. optimize |

1526 | integer iwm |

1527 | integer iownd, iowns, |

1528 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1529 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1530 | integer i, meband, ml, mu |

1531 | double precision wm, x, tem |

1532 | double precision rowns, |

1533 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |

1534 | double precision di, hl0, phl0, r |

1535 | dimension wm(1), iwm(1), x(1), tem(1) |

1536 | common /ls0001/ rowns(209), |

1537 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |

1538 | 3 iownd(14), iowns(6), |

1539 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1540 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1541 | c----------------------------------------------------------------------- |

1542 | c this routine manages the solution of the linear system arising from |

1543 | c a chord iteration. it is called if miter .ne. 0. |

1544 | c if miter is 1 or 2, it calls dgesl to accomplish this. |

1545 | c if miter = 3 it updates the coefficient h*el0 in the diagonal |

1546 | c matrix, and then computes the solution. |

1547 | c if miter is 4 or 5, it calls dgbsl. |

1548 | c communication with solsy uses the following variables.. |

1549 | c wm = real work space containing the inverse diagonal matrix if |

1550 | c miter = 3 and the lu decomposition of the matrix otherwise. |

1551 | c storage of matrix elements starts at wm(3). |

1552 | c wm also contains the following matrix-related data.. |

1553 | c wm(1) = sqrt(uround) (not used here), |

1554 | c wm(2) = hl0, the previous value of h*el0, used if miter = 3. |

1555 | c iwm = integer work space containing pivot information, starting at |

1556 | c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band |

1557 | c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. |

1558 | c x = the right-hand side vector on input, and the solution vector |

1559 | c on output, of length n. |

1560 | c tem = vector of work space of length n, not used in this version. |

1561 | c iersl = output flag (in common). iersl = 0 if no trouble occurred. |

1562 | c iersl = 1 if a singular matrix arose with miter = 3. |

1563 | c this routine also uses the common variables el0, h, miter, and n. |

1564 | c----------------------------------------------------------------------- |

1565 | iersl = 0 |

1566 | go to (100, 100, 300, 400, 400), miter |

1567 | 100 call dgesl (wm(3), n, n, iwm(21), x, 0) |

1568 | return |

1569 | c |

1570 | 300 phl0 = wm(2) |

1571 | hl0 = h*el0 |

1572 | wm(2) = hl0 |

1573 | if (hl0 .eq. phl0) go to 330 |

1574 | r = hl0/phl0 |

1575 | do 320 i = 1,n |

1576 | di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2)) |

1577 | if (dabs(di) .eq. 0.0d0) go to 390 |

1578 | 320 wm(i+2) = 1.0d0/di |

1579 | 330 do 340 i = 1,n |

1580 | 340 x(i) = wm(i+2)*x(i) |

1581 | return |

1582 | 390 iersl = 1 |

1583 | return |

1584 | c |

1585 | 400 ml = iwm(1) |

1586 | mu = iwm(2) |

1587 | meband = 2*ml + mu + 1 |

1588 | call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0) |

1589 | return |

1590 | c----------------------- end of subroutine solsy ----------------------- |

1591 | end |

1592 | subroutine prepj (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm, |

1593 | 1 f, jac) |

1594 | clll. optimize |

1595 | cascend changes |

1596 | external f, jac |

1597 | integer neq, nyh, iwm |

1598 | integer iownd, iowns, |

1599 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1600 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1601 | integer i, i1, i2, ier, ii, j, j1, jj, lenp, |

1602 | 1 mba, mband, meb1, meband, ml, ml3, mu, np1 |

1603 | double precision y, yh, ewt, ftem, savf, wm |

1604 | double precision rowns, |

1605 | 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |

1606 | double precision con, di, fac, hl0, r, r0, srur, yi, yj, yjj, |

1607 | 1 vnorm |

1608 | Ckaa double precision time1, time2 |

1609 | dimension neq(1), y(1), yh(nyh,1), ewt(1), ftem(1), savf(1), |

1610 | 1 wm(1), iwm(1) |

1611 | common /ls0001/ rowns(209), |

1612 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, |

1613 | 3 iownd(14), iowns(6), |

1614 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1615 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1616 | c----------------------------------------------------------------------- |

1617 | c prepj is called by stode to compute and process the matrix |

1618 | c p = i - h*el(1)*j , where j is an approximation to the jacobian. |

1619 | c here j is computed by the user-supplied routine jac if |

1620 | c miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5. |

1621 | c if miter = 3, a diagonal approximation to j is used. |

1622 | c j is stored in wm and replaced by p. if miter .ne. 3, p is then |

1623 | c subjected to lu decomposition in preparation for later solution |

1624 | c of linear systems with p as coefficient matrix. this is done |

1625 | c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5. |

1626 | c |

1627 | c in addition to variables described previously, communication |

1628 | c with prepj uses the following.. |

1629 | c y = array containing predicted values on entry. |

1630 | c ftem = work array of length n (acor in stode). |

1631 | c savf = array containing f evaluated at predicted y. |

1632 | c wm = real work space for matrices. on output it contains the |

1633 | c inverse diagonal matrix if miter = 3 and the lu decomposition |

1634 | c of p if miter is 1, 2 , 4, or 5. |

1635 | c storage of matrix elements starts at wm(3). |

1636 | c wm also contains the following matrix-related data.. |

1637 | c wm(1) = sqrt(uround), used in numerical jacobian increments. |

1638 | c wm(2) = h*el0, saved for later use if miter = 3. |

1639 | c iwm = integer work space containing pivot information, starting at |

1640 | c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band |

1641 | c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. |

1642 | c el0 = el(1) (input). |

1643 | c ierpj = output error flag, = 0 if no trouble, .gt. 0 if |

1644 | c p matrix found to be singular. |

1645 | c jcur = output flag = 1 to indicate that the jacobian matrix |

1646 | c (or approximation) is now current. |

1647 | c this routine also uses the common variables el0, h, tn, uround, |

1648 | c miter, n, nfe, and nje. |

1649 | c----------------------------------------------------------------------- |

1650 | nje = nje + 1 |

1651 | ierpj = 0 |

1652 | jcur = 1 |

1653 | hl0 = h*el0 |

1654 | go to (100, 200, 300, 400, 500), miter |

1655 | c if miter = 1, call jac and multiply by scalar. ----------------------- |

1656 | 100 lenp = n*n |

1657 | do 110 i = 1,lenp |

1658 | 110 wm(i+2) = 0.0d0 |

1659 | call jac (neq, tn, y, 0, 0, wm(3), n) |

1660 | con = -hl0 |

1661 | do 120 i = 1,lenp |

1662 | 120 wm(i+2) = wm(i+2)*con |

1663 | go to 240 |

1664 | c if miter = 2, make n calls to f to approximate j. -------------------- |

1665 | 200 fac = vnorm (n, savf, ewt) |

1666 | r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac |

1667 | if (r0 .eq. 0.0d0) r0 = 1.0d0 |

1668 | srur = wm(1) |

1669 | j1 = 2 |

1670 | do 230 j = 1,n |

1671 | yj = y(j) |

1672 | r = dmax1(srur*dabs(yj),r0/ewt(j)) |

1673 | y(j) = y(j) + r |

1674 | fac = -hl0/r |

1675 | call f (neq, tn, y, ftem) |

1676 | do 220 i = 1,n |

1677 | 220 wm(i+j1) = (ftem(i) - savf(i))*fac |

1678 | y(j) = yj |

1679 | j1 = j1 + n |

1680 | 230 continue |

1681 | nfe = nfe + n |

1682 | c add identity matrix. ------------------------------------------------- |

1683 | 240 j = 3 |

1684 | np1 = n + 1 |

1685 | do 250 i = 1,n |

1686 | wm(j) = wm(j) + 1.0d0 |

1687 | 250 j = j + np1 |

1688 | c do lu decomposition on p. -------------------------------------------- |

1689 | Ckaa call aftime(time1) |

1690 | call dgefa (wm(3), n, n, iwm(21), ier) |

1691 | Ckaa call aftime(time2) |

1692 | Ckaa time2 = time2 - time1 |

1693 | Ckaa write (6,*) "Time for decomposition" |

1694 | Ckaa write (6,700) time2 |

1695 | if (ier .ne. 0) ierpj = 1 |

1696 | return |

1697 | c if miter = 3, construct a diagonal approximation to j and p. --------- |

1698 | 300 wm(2) = hl0 |

1699 | r = el0*0.1d0 |

1700 | do 310 i = 1,n |

1701 | 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2)) |

1702 | call f (neq, tn, y, wm(3)) |

1703 | nfe = nfe + 1 |

1704 | do 320 i = 1,n |

1705 | r0 = h*savf(i) - yh(i,2) |

1706 | di = 0.1d0*r0 - h*(wm(i+2) - savf(i)) |

1707 | wm(i+2) = 1.0d0 |

1708 | if (dabs(r0) .lt. uround/ewt(i)) go to 320 |

1709 | if (dabs(di) .eq. 0.0d0) go to 330 |

1710 | wm(i+2) = 0.1d0*r0/di |

1711 | 320 continue |

1712 | return |

1713 | 330 ierpj = 1 |

1714 | return |

1715 | c if miter = 4, call jac and multiply by scalar. ----------------------- |

1716 | 400 ml = iwm(1) |

1717 | mu = iwm(2) |

1718 | ml3 = ml + 3 |

1719 | mband = ml + mu + 1 |

1720 | meband = mband + ml |

1721 | lenp = meband*n |

1722 | do 410 i = 1,lenp |

1723 | 410 wm(i+2) = 0.0d0 |

1724 | call jac (neq, tn, y, ml, mu, wm(ml3), meband) |

1725 | con = -hl0 |

1726 | do 420 i = 1,lenp |

1727 | 420 wm(i+2) = wm(i+2)*con |

1728 | go to 570 |

1729 | c if miter = 5, make mband calls to f to approximate j. ---------------- |

1730 | 500 ml = iwm(1) |

1731 | mu = iwm(2) |

1732 | mband = ml + mu + 1 |

1733 | mba = min0(mband,n) |

1734 | meband = mband + ml |

1735 | meb1 = meband - 1 |

1736 | srur = wm(1) |

1737 | fac = vnorm (n, savf, ewt) |

1738 | r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac |

1739 | if (r0 .eq. 0.0d0) r0 = 1.0d0 |

1740 | do 560 j = 1,mba |

1741 | do 530 i = j,n,mband |

1742 | yi = y(i) |

1743 | r = dmax1(srur*dabs(yi),r0/ewt(i)) |

1744 | 530 y(i) = y(i) + r |

1745 | call f (neq, tn, y, ftem) |

1746 | do 550 jj = j,n,mband |

1747 | y(jj) = yh(jj,1) |

1748 | yjj = y(jj) |

1749 | r = dmax1(srur*dabs(yjj),r0/ewt(jj)) |

1750 | fac = -hl0/r |

1751 | i1 = max0(jj-mu,1) |

1752 | i2 = min0(jj+ml,n) |

1753 | ii = jj*meb1 - ml + 2 |

1754 | do 540 i = i1,i2 |

1755 | 540 wm(ii+i) = (ftem(i) - savf(i))*fac |

1756 | 550 continue |

1757 | 560 continue |

1758 | nfe = nfe + mba |

1759 | c add identity matrix. ------------------------------------------------- |

1760 | 570 ii = mband + 2 |

1761 | do 580 i = 1,n |

1762 | wm(ii) = wm(ii) + 1.0d0 |

1763 | 580 ii = ii + meband |

1764 | c do lu decomposition of p. -------------------------------------------- |

1765 | call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier) |

1766 | if (ier .ne. 0) ierpj = 1 |

1767 | return |

1768 | 700 format(d21.13) |

1769 | c----------------------- end of subroutine prepj ----------------------- |

1770 | end |

1771 | subroutine stode (neq, y, yh, nyh, yh1, ewt, savf, acor, |

1772 | 1 wm, iwm, f, jac, pjac, slvs) |

1773 | clll. optimize |

1774 | external f, jac, pjac, slvs |

1775 | integer neq, nyh, iwm |

1776 | integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp, |

1777 | 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1778 | 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1779 | integer i, i1, iredo, iret, j, jb, m, ncf, newq |

1780 | double precision y, yh, yh1, ewt, savf, acor, wm |

1781 | double precision conit, crate, el, elco, hold, rmax, tesco, |

1782 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround |

1783 | double precision dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, |

1784 | 1 r, rh, rhdn, rhsm, rhup, told, vnorm |

1785 | dimension neq(1), y(1), yh(nyh,1), yh1(1), ewt(1), savf(1), |

1786 | 1 acor(1), wm(1), iwm(1) |

1787 | common /ls0001/ conit, crate, el(13), elco(13,12), |

1788 | 1 hold, rmax, tesco(3,12), |

1789 | 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14), |

1790 | 3 ialth, ipup, lmax, meo, nqnyh, nslp, |

1791 | 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, |

1792 | 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu |

1793 | c----------------------------------------------------------------------- |

1794 | c stode performs one step of the integration of an initial value |

1795 | c problem for a system of ordinary differential equations. |

1796 | c note.. stode is independent of the value of the iteration method |

1797 | c indicator miter, when this is .ne. 0, and hence is independent |

1798 | c of the type of chord method used, or the jacobian structure. |

1799 | c communication with stode is done with the following variables.. |

1800 | c |

1801 | c neq = integer array containing problem size in neq(1), and |

1802 | c passed as the neq argument in all calls to f and jac. |

1803 | c y = an array of length .ge. n used as the y argument in |

1804 | c all calls to f and jac. |

1805 | c yh = an nyh by lmax array containing the dependent variables |

1806 | c and their approximate scaled derivatives, where |

1807 | c lmax = maxord + 1. yh(i,j+1) contains the approximate |

1808 | c j-th derivative of y(i), scaled by h**j/factorial(j) |

1809 | c (j = 0,1,...,nq). on entry for the first step, the first |

1810 | c two columns of yh must be set from the initial values. |

1811 | c nyh = a constant integer .ge. n, the first dimension of yh. |

1812 | c yh1 = a one-dimensional array occupying the same space as yh. |

1813 | c ewt = an array of length n containing multiplicative weights |

1814 | c for local error measurements. local errors in y(i) are |

1815 | c compared to 1.0/ewt(i) in various error tests. |

1816 | c savf = an array of working storage, of length n. |

1817 | c also used for input of yh(*,maxord+2) when jstart = -1 |

1818 | c and maxord .lt. the current order nq. |

1819 | c acor = a work array of length n, used for the accumulated |

1820 | c corrections. on a successful return, acor(i) contains |

1821 | c the estimated one-step local error in y(i). |

1822 | c wm,iwm = real and integer work arrays associated with matrix |

1823 | c operations in chord iteration (miter .ne. 0). |

1824 | c pjac = name of routine to evaluate and preprocess jacobian matrix |

1825 | c and p = i - h*el0*jac, if a chord method is being used. |

1826 | c slvs = name of routine to solve linear system in chord iteration. |

1827 | c ccmax = maximum relative change in h*el0 before pjac is called. |

1828 | c h = the step size to be attempted on the next step. |

1829 | c h is altered by the error control algorithm during the |

1830 | c problem. h can be either positive or negative, but its |

1831 | c sign must remain constant throughout the problem. |

1832 | c hmin = the minimum absolute value of the step size h to be used. |

1833 | c hmxi = inverse of the maximum absolute value of h to be used. |

1834 | c hmxi = 0.0 is allowed and corresponds to an infinite hmax. |

1835 | c hmin and hmxi may be changed at any time, but will not |

1836 | c take effect until the next change of h is considered. |

1837 | c tn = the independent variable. tn is updated on each step taken. |

1838 | c jstart = an integer used for input only, with the following |

1839 | c values and meanings.. |

1840 | c 0 perform the first step. |

1841 | c .gt.0 take a new step continuing from the last. |

1842 | c -1 take the next step with a new value of h, maxord, |

1843 | c n, meth, miter, and/or matrix parameters. |

1844 | c -2 take the next step with a new value of h, |

1845 | c but with other inputs unchanged. |

1846 | c on return, jstart is set to 1 to facilitate continuation. |

1847 | c kflag = a completion code with the following meanings.. |

1848 | c 0 the step was succesful. |

1849 | c -1 the requested error could not be achieved. |

1850 | c -2 corrector convergence could not be achieved. |

1851 | c -3 fatal error in pjac or slvs. |

1852 | c a return with kflag = -1 or -2 means either |

1853 | c abs(h) = hmin or 10 consecutive failures occurred. |

1854 | c on a return with kflag negative, the values of tn and |

1855 | c the yh array are as of the beginning of the last |

1856 | c step, and h is the last step size attempted. |

1857 | c maxord = the maximum order of integration method to be allowed. |

1858 | c maxcor = the maximum number of corrector iterations allowed. |

1859 | c msbp = maximum number of steps between pjac calls (miter .gt. 0). |

1860 | c mxncf = maximum number of convergence failures allowed. |

1861 | c meth/miter = the method flags. see description in driver. |

1862 | c n = the number of first-order differential equations. |

1863 | c----------------------------------------------------------------------- |

1864 | kflag = 0 |

1865 | told = tn |

1866 | ncf = 0 |

1867 | ierpj = 0 |

1868 | iersl = 0 |

1869 | jcur = 0 |

1870 | icf = 0 |

1871 | delp = 0.0d0 |

1872 | if (jstart .gt. 0) go to 200 |

1873 | if (jstart .eq. -1) go to 100 |

1874 | if (jstart .eq. -2) go to 160 |

1875 | c----------------------------------------------------------------------- |

1876 | c on the first call, the order is set to 1, and other variables are |

1877 | c initialized. rmax is the maximum ratio by which h can be increased |

1878 | c in a single step. it is initially 1.e4 to compensate for the small |

1879 | c initial h, but then is normally equal to 10. if a failure |

1880 | c occurs (in corrector convergence or error test), rmax is set at 2 |

1881 | c for the next increase. |

1882 | c----------------------------------------------------------------------- |

1883 | lmax = maxord + 1 |

1884 | nq = 1 |

1885 | l = 2 |

1886 | ialth = 2 |

1887 | rmax = 10000.0d0 |

1888 | rc = 0.0d0 |

1889 | el0 = 1.0d0 |

1890 | crate = 0.7d0 |

1891 | hold = h |

1892 | meo = meth |

1893 | nslp = 0 |

1894 | ipup = miter |

1895 | iret = 3 |

1896 | go to 140 |

1897 | c----------------------------------------------------------------------- |

1898 | c the following block handles preliminaries needed when jstart = -1. |

1899 | c ipup is set to miter to force a matrix update. |

1900 | c if an order increase is about to be considered (ialth = 1), |

1901 | c ialth is reset to 2 to postpone consideration one more step. |

1902 | c if the caller has changed meth, cfode is called to reset |

1903 | c the coefficients of the method. |

1904 | c if the caller has changed maxord to a value less than the current |

1905 | c order nq, nq is reduced to maxord, and a new h chosen accordingly. |

1906 | c if h is to be changed, yh must be rescaled. |

1907 | c if h or meth is being changed, ialth is reset to l = nq + 1 |

1908 | c to prevent further changes in h for that many steps. |

1909 | c----------------------------------------------------------------------- |

1910 | 100 ipup = miter |

1911 | lmax = maxord + 1 |

1912 | if (ialth .eq. 1) ialth = 2 |

1913 | if (meth .eq. meo) go to 110 |

1914 | call cfode (meth, elco, tesco) |

1915 | meo = meth |

1916 | if (nq .gt. maxord) go to 120 |

1917 | ialth = l |

1918 | iret = 1 |

1919 | go to 150 |

1920 | 110 if (nq .le. maxord) go to 160 |

1921 | 120 nq = maxord |

1922 | l = lmax |

1923 | do 125 i = 1,l |

1924 | 125 el(i) = elco(i,nq) |

1925 | nqnyh = nq*nyh |

1926 | rc = rc*el(1)/el0 |

1927 | el0 = el(1) |

1928 | conit = 0.5d0/dfloat(nq+2) |

1929 | ddn = vnorm (n, savf, ewt)/tesco(1,l) |

1930 | exdn = 1.0d0/dfloat(l) |

1931 | rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0) |

1932 | rh = dmin1(rhdn,1.0d0) |

1933 | iredo = 3 |

1934 | if (h .eq. hold) go to 170 |

1935 | rh = dmin1(rh,dabs(h/hold)) |

1936 | h = hold |

1937 | go to 175 |

1938 | c----------------------------------------------------------------------- |

1939 | c cfode is called to get all the integration coefficients for the |

1940 | c current meth. then the el vector and related constants are reset |

1941 | c whenever the order nq is changed, or at the start of the problem. |

1942 | c----------------------------------------------------------------------- |

1943 | 140 call cfode (meth, elco, tesco) |

1944 | 150 do 155 i = 1,l |

1945 | 155 el(i) = elco(i,nq) |

1946 |