Parent Directory | Revision Log

Revision **1** -
(**hide annotations**)
(**download**)

*Fri Oct 29 20:54:12 2004 UTC*
(19 years, 4 months ago)
by *aw0a*

Original Path:*trunk/lsod/lsode.f*

File size: 118049 byte(s)

Original Path:

File size: 118049 byte(s)

Setting up web subdirectory in repository

1 | aw0a | 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 = |