1 subroutine lsode (f, neq, y, t, tout, itol, rtol, atol, itask, 2 1 istate, iopt, rwork, lrw, iwork, liw, jac, mf) 3 4 external f, jac 5 integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf 6 double precision y, t, tout, rtol, atol, rwork 7 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw) 8c 9c!purpose 10c livermore solver for ordinary differential equations. 11c this version is in double precision. 12c 13c lsode solves the initial value problem for stiff or nonstiff 14c systems of first order ode-s, 15c dy/dt = f(t,y) , or, in component form, 16c dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq). 17c lsode is a package based on the gear and gearb packages, and on the 18c october 23, 1978 version of the tentative odepack user interface 19c standard, with minor modifications. 20c 21c!summary of usage. 22c 23c communication between the user and the lsode package, for normal 24c situations, is summarized here. this summary describes only a subset 25c of the full set of options available. see the full description for 26c details, including optional communication, nonstandard options, 27c and instructions for special situations. see also the example 28c problem (with program and output) following this summary. 29c 30c a. first provide a subroutine of the form.. 31c subroutine f (neq, t, y, ydot) 32c dimension y(neq), ydot(neq) 33c which supplies the vector function f by loading ydot(i) with f(i). 34c 35c b. next determine (or guess) whether or not the problem is stiff. 36c stiffness occurs when the jacobian matrix df/dy has an eigenvalue 37c whose real part is negative and large in magnitude, compared to the 38c reciprocal of the t span of interest. if the problem is nonstiff, 39c use a method flag mf = 10. if it is stiff, there are four standard 40c choices for mf, and lsode requires the jacobian matrix in some form. 41c this matrix is regarded either as full (mf = 21 or 22), 42c or banded (mf = 24 or 25). in the banded case, lsode requires two 43c half-bandwidth parameters ml and mu. these are, respectively, the 44c widths of the lower and upper parts of the band, excluding the main 45c diagonal. thus the band consists of the locations (i,j) with 46c i-ml .le. j .le. i+mu, and the full bandwidth is ml+mu+1. 47c 48c c. if the problem is stiff, you are encouraged to supply the jacobian 49c directly (mf = 21 or 24), but if this is not feasible, lsode will 50c compute it internally by difference quotients (mf = 22 or 25). 51c if you are supplying the jacobian, provide a subroutine of the form.. 52c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) 53c dimension y(neq), pd(nrowpd,neq) 54c which supplies df/dy by loading pd as follows.. 55c for a full jacobian (mf = 21), load pd(i,j) with df(i)/dy(j), 56c the partial derivative of f(i) with respect to y(j). (ignore the 57c ml and mu arguments in this case.) 58c for a banded jacobian (mf = 24), load pd(i-j+mu+1,j) with 59c df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of 60c pd from the top down. 61c in either case, only nonzero elements need be loaded. 62c 63c d. write a main program which calls subroutine lsode once for 64c each point at which answers are desired. this should also provide 65c for possible use of logical unit 6 for output of error messages 66c by lsode. on the first call to lsode, supply arguments as follows.. 67c f = name of subroutine for right-hand side vector f. 68c this name must be declared external in calling program. 69c neq = number of first order ode-s. 70c y = array of initial values, of length neq. 71c t = the initial value of the independent variable. 72c tout = first point where output is desired (.ne. t). 73c itol = 1 or 2 according as atol (below) is a scalar or array. 74c rtol = relative tolerance parameter (scalar). 75c atol = absolute tolerance parameter (scalar or array). 76c the estimated local error in y(i) will be controlled so as 77c to be roughly less (in magnitude) than 78c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or 79c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2. 80c thus the local error test passes if, in each component, 81c either the absolute error is less than atol (or atol(i)), 82c or the relative error is less than rtol. 83c use rtol = 0.0 for pure absolute error control, and 84c use atol = 0.0 (or atol(i) = 0.0) for pure relative error 85c control. caution.. actual (global) errors may exceed these 86c local tolerances, so choose them conservatively. 87c itask = 1 for normal computation of output values of y at t = tout. 88c istate = integer flag (input and output). set istate = 1. 89c iopt = 0 to indicate no optional inputs used. 90c rwork = real work array of length at least.. 91c 20 + 16*neq for mf = 10, 92c 22 + 9*neq + neq**2 for mf = 21 or 22, 93c 22 + 10*neq + (2*ml + mu)*neq for mf = 24 or 25. 94c lrw = declared length of rwork (in user-s dimension). 95c iwork = integer work array of length at least.. 96c 20 for mf = 10, 97c 20 + neq for mf = 21, 22, 24, or 25. 98c if mf = 24 or 25, input in iwork(1),iwork(2) the lower 99c and upper half-bandwidths ml,mu. 100c liw = declared length of iwork (in user-s dimension). 101c jac = name of subroutine for jacobian matrix (mf = 21 or 24). 102c if used, this name must be declared external in calling 103c program. if not used, pass a dummy name. 104c mf = method flag. standard values are.. 105c 10 for nonstiff (adams) method, no jacobian used. 106c 21 for stiff (bdf) method, user-supplied full jacobian. 107c 22 for stiff method, internally generated full jacobian. 108c 24 for stiff method, user-supplied banded jacobian. 109c 25 for stiff method, internally generated banded jacobian. 110c note that the main program must declare arrays y, rwork, iwork, 111c and possibly atol. 112c 113c e. the output from the first call (or any call) is.. 114c y = array of computed values of y(t) vector. 115c t = corresponding value of independent variable (normally tout). 116c istate = 2 if lsode was successful, negative otherwise. 117c -1 means excess work done on this call (perhaps wrong mf). 118c -2 means excess accuracy requested (tolerances too small). 119c -3 means illegal input detected (see printed message). 120c -4 means repeated error test failures (check all inputs). 121c -5 means repeated convergence failures (perhaps bad jacobian 122c supplied or wrong choice of mf or tolerances). 123c -6 means error weight became zero during problem. (solution 124c component i vanished, and atol or atol(i) = 0.) 125c 126c f. to continue the integration after a successful return, simply 127c reset tout and call lsode again. no other parameters need be reset. 128c 129c 130c!example problem. 131c 132c the following is a simple example problem, with the coding 133c needed for its solution by lsode. the problem is from chemical 134c kinetics, and consists of the following three rate equations.. 135c dy1/dt = -.04*y1 + 1.e4*y2*y3 136c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2 137c dy3/dt = 3.e7*y2**2 138c on the interval from t = 0.0 to t = 4.e10, with initial conditions 139c y1 = 1.0, y2 = y3 = 0. the problem is stiff. 140c 141c the following coding solves this problem with lsode, using mf = 21 142c and printing results at t = .4, 4., ..., 4.e10. it uses 143c itol = 2 and atol much smaller for y2 than y1 or y3 because 144c y2 has much smaller values. 145c at the end of the run, statistical quantities of interest are 146c printed (see optional outputs in the full description below). 147c 148c external fex, jex 149c double precision atol, rwork, rtol, t, tout, y 150c dimension y(3), atol(3), rwork(58), iwork(23) 151c neq = 3 152c y(1) = 1.0d+0 153c y(2) = 0.0d+0 154c y(3) = 0.0d+0 155c t = 0.0d+0 156c tout = .40d+0 157c itol = 2 158c rtol = 1.0d-4 159c atol(1) = 1.0d-6 160c atol(2) = 1.0d-10 161c atol(3) = 1.0d-6 162c itask = 1 163c istate = 1 164c iopt = 0 165c lrw = 58 166c liw = 23 167c mf = 21 168c do 40 iout = 1,12 169c call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate, 170c 1 iopt,rwork,lrw,iwork,liw,jex,mf) 171c write(6,20)t,y(1),y(2),y(3) 172c 20 format(7h at t =,e12.4,6h y =,3e14.6) 173c if (istate .lt. 0) go to 80 174c 40 tout = tout*10.0d+0 175c write(6,60)iwork(11),iwork(12),iwork(13) 176c 60 format(/12h no. steps =,i4,11h no. f-s =,i4,11h no. j-s =,i4) 177c stop 178c 80 write(6,90)istate 179c 90 format(///22h error halt.. istate =,i3) 180c stop 181c end 182c 183c subroutine fex (neq, t, y, ydot) 184c double precision t, y, ydot 185c dimension y(3), ydot(3) 186c ydot(1) = -.040d+0*y(1) + 1.0d+4*y(2)*y(3) 187c ydot(3) = 3.0d+7*y(2)*y(2) 188c ydot(2) = -ydot(1) - ydot(3) 189c return 190c end 191c 192c subroutine jex (neq, t, y, ml, mu, pd, nrpd) 193c double precision pd, t, y 194c dimension y(3), pd(nrpd,3) 195c pd(1,1) = -0.040d+0 196c pd(1,2) = 1.0d+4*y(3) 197c pd(1,3) = 1.0d+4*y(2) 198c pd(2,1) = 0.040d+0 199c pd(2,3) = -pd(1,3) 200c pd(3,2) = 6.0d+7*y(2) 201c pd(2,2) = -pd(1,2) - pd(3,2) 202c return 203c end 204c 205c the output of this program (on a cdc-7600 in single precision) 206c is as follows.. 207c 208c at t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02 209c at t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02 210c at t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01 211c at t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01 212c at t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01 213c at t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01 214c at t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01 215c at t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01 216c at t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01 217c at t = 4.0000e+08 y = 5.494529e-06 2.197824e-11 9.999945e-01 218c at t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01 219c at t = 4.0000e+10 y = -7.170586e-08 -2.868234e-13 1.000000e+00 220c 221c no. steps = 330 no. f-s = 405 no. j-s = 69 222c 223c!full description of user interface to lsode. 224c 225c the user interface to lsode consists of the following parts. 226c 227c i. the call sequence to subroutine lsode, which is a driver 228c routine for the solver. this includes descriptions of both 229c the call sequence arguments and of user-supplied routines. 230c following these descriptions is a description of 231c optional inputs available through the call sequence, and then 232c a description of optional outputs (in the work arrays). 233c 234c ii. descriptions of other routines in the lsode package that may be 235c (optionally) called by the user. these provide the ability to 236c alter error message handling, save and restore the internal 237c common, and obtain specified derivatives of the solution y(t). 238c 239c iii. descriptions of common blocks to be declared in overlay 240c or similar environments, or to be saved when doing an interrupt 241c of the problem and continued solution later. 242c 243c iv. description of two subroutines in the lsode package, either of 244c which the user may replace with his own version, if desired. 245c these relate to the measurement of errors. 246c 247c 248c part i. call sequence. 249c 250c the call sequence parameters used for input only are 251c f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf, 252c and those used for both input and output are 253c y, t, istate. 254c the work arrays rwork and iwork are also used for conditional and 255c optional inputs and optional outputs. (the term output here refers 256c to the return from subroutine lsode to the user-s calling program.) 257c 258c the legality of input parameters will be thoroughly checked on the 259c initial call for the problem, but not checked thereafter unless a 260c change in input parameters is flagged by istate = 3 on input. 261c 262c the descriptions of the call arguments are as follows. 263c 264c f = the name of the user-supplied subroutine defining the 265c ode system. the system must be put in the first-order 266c form dy/dt = f(t,y), where f is a vector-valued function 267c of the scalar t and the vector y. subroutine f is to 268c compute the function f. it is to have the form 269c subroutine f (neq, t, y, ydot) 270c dimension y(1), ydot(1) 271c where neq, t, and y are input, and the array ydot = f(t,y) 272c is output. y and ydot are arrays of length neq. 273c (in the dimension statement above, 1 is a dummy 274c dimension.. it can be replaced by any value.) 275c subroutine f should not alter y(1),...,y(neq). 276c f must be declared external in the calling program. 277c 278c subroutine f may access user-defined quantities in 279c neq(2),... and y(neq(1)+1),... if neq is an array 280c (dimensioned in f) and y has length exceeding neq(1). 281c see the descriptions of neq and y below. 282c 283c neq = the size of the ode system (number of first order 284c ordinary differential equations). used only for input. 285c neq may be decreased, but not increased, during the problem. 286c if neq is decreased (with istate = 3 on input), the 287c remaining components of y should be left undisturbed, if 288c these are to be accessed in f and/or jac. 289c 290c normally, neq is a scalar, and it is generally referred to 291c as a scalar in this user interface description. however, 292c neq may be an array, with neq(1) set to the system size. 293c (the lsode package accesses only neq(1).) in either case, 294c this parameter is passed as the neq argument in all calls 295c to f and jac. hence, if it is an array, locations 296c neq(2),... may be used to store other integer data and pass 297c it to f and/or jac. subroutines f and/or jac must include 298c neq in a dimension statement in that case. 299c 300c y = a real array for the vector of dependent variables, of 301c length neq or more. used for both input and output on the 302c first call (istate = 1), and only for output on other calls. 303c on the first call, y must contain the vector of initial 304c values. on output, y contains the computed solution vector, 305c evaluated at t. if desired, the y array may be used 306c for other purposes between calls to the solver. 307c 308c this array is passed as the y argument in all calls to 309c f and jac. hence its length may exceed neq, and locations 310c y(neq+1),... may be used to store other real data and 311c pass it to f and/or jac. (the lsode package accesses only 312c y(1),...,y(neq).) 313c 314c t = the independent variable. on input, t is used only on the 315c first call, as the initial point of the integration. 316c on output, after each call, t is the value at which a 317c computed solution y is evaluated (usually the same as tout). 318c on an error return, t is the farthest point reached. 319c 320c tout = the next value of t at which a computed solution is desired. 321c used only for input. 322c 323c when starting the problem (istate = 1), tout may be equal 324c to t for one call, then should .ne. t for the next call. 325c for the initial t, an input value of tout .ne. t is used 326c in order to determine the direction of the integration 327c (i.e. the algebraic sign of the step sizes) and the rough 328c scale of the problem. integration in either direction 329c (forward or backward in t) is permitted. 330c 331c if itask = 2 or 5 (one-step modes), tout is ignored after 332c the first call (i.e. the first call with tout .ne. t). 333c otherwise, tout is required on every call. 334c 335c if itask = 1, 3, or 4, the values of tout need not be 336c monotone, but a value of tout which backs up is limited 337c to the current internal t interval, whose endpoints are 338c tcur - hu and tcur (see optional outputs, below, for 339c tcur and hu). 340c 341c itol = an indicator for the type of error control. see 342c description below under atol. used only for input. 343c 344c rtol = a relative error tolerance parameter, either a scalar or 345c an array of length neq. see description below under atol. 346c input only. 347c 348c atol = an absolute error tolerance parameter, either a scalar or 349c an array of length neq. input only. 350c 351c the input parameters itol, rtol, and atol determine 352c the error control performed by the solver. the solver will 353c control the vector e = (e(i)) of estimated local errors 354c in y, according to an inequality of the form 355c rms-norm of ( e(i)/ewt(i) ) .le. 1, 356c where ewt(i) = rtol(i)*abs(y(i)) + atol(i), 357c and the rms-norm (root-mean-square norm) here is 358c rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i)) 359c is a vector of weights which must always be positive, and 360c the values of rtol and atol should all be non-negative. 361c the following table gives the types (scalar/array) of 362c rtol and atol, and the corresponding form of ewt(i). 363c 364c itol rtol atol ewt(i) 365c 1 scalar scalar rtol*abs(y(i)) + atol 366c 2 scalar array rtol*abs(y(i)) + atol(i) 367c 3 array scalar rtol(i)*abs(y(i)) + atol 368c 4 array array rtol(i)*abs(y(i)) + atol(i) 369c 370c when either of these parameters is a scalar, it need not 371c be dimensioned in the user-s calling program. 372c 373c if none of the above choices (with itol, rtol, and atol 374c fixed throughout the problem) is suitable, more general 375c error controls can be obtained by substituting 376c user-supplied routines for the setting of ewt and/or for 377c the norm calculation. see part iv below. 378c 379c if global errors are to be estimated by making a repeated 380c run on the same problem with smaller tolerances, then all 381c components of rtol and atol (i.e. of ewt) should be scaled 382c down uniformly. 383c 384c itask = an index specifying the task to be performed. 385c input only. itask has the following values and meanings. 386c 1 means normal computation of output values of y(t) at 387c t = tout (by overshooting and interpolating). 388c 2 means take one step only and return. 389c 3 means stop at the first internal mesh point at or 390c beyond t = tout and return. 391c 4 means normal computation of output values of y(t) at 392c t = tout but without overshooting t = tcrit. 393c tcrit must be input as rwork(1). tcrit may be equal to 394c or beyond tout, but not behind it in the direction of 395c integration. this option is useful if the problem 396c has a singularity at or beyond t = tcrit. 397c 5 means take one step, without passing tcrit, and return. 398c tcrit must be input as rwork(1). 399c 400c note.. if itask = 4 or 5 and the solver reaches tcrit 401c (within roundoff), it will return t = tcrit (exactly) to 402c indicate this (unless itask = 4 and tout comes before tcrit, 403c in which case answers at t = tout are returned first). 404c 405c istate = an index used for input and output to specify the 406c the state of the calculation. 407c 408c on input, the values of istate are as follows. 409c 1 means this is the first call for the problem 410c (initializations will be done). see note below. 411c 2 means this is not the first call, and the calculation 412c is to continue normally, with no change in any input 413c parameters except possibly tout and itask. 414c (if itol, rtol, and/or atol are changed between calls 415c with istate = 2, the new values will be used but not 416c tested for legality.) 417c 3 means this is not the first call, and the 418c calculation is to continue normally, but with 419c a change in input parameters other than 420c tout and itask. changes are allowed in 421c neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu, 422c and any of the optional inputs except h0. 423c (see iwork description for ml and mu.) 424c note.. a preliminary call with tout = t is not counted 425c as a first call here, as no initialization or checking of 426c input is done. (such a call is sometimes useful for the 427c purpose of outputting the initial conditions.) 428c thus the first call for which tout .ne. t requires 429c istate = 1 on input. 430c 431c on output, istate has the following values and meanings. 432c 1 means nothing was done, as tout was equal to t with 433c istate = 1 on input. (however, an internal counter was 434c set to detect and prevent repeated calls of this type.) 435c 2 means the integration was performed successfully. 436c -1 means an excessive amount of work (more than mxstep 437c steps) was done on this call, before completing the 438c requested task, but the integration was otherwise 439c successful as far as t. (mxstep is an optional input 440c and is normally 500.) to continue, the user may 441c simply reset istate to a value .gt. 1 and call again 442c (the excess work step counter will be reset to 0). 443c in addition, the user may increase mxstep to avoid 444c this error return (see below on optional inputs). 445c -2 means too much accuracy was requested for the precision 446c of the machine being used. this was detected before 447c completing the requested task, but the integration 448c was successful as far as t. to continue, the tolerance 449c parameters must be reset, and istate must be set 450c to 3. the optional output tolsf may be used for this 451c purpose. (note.. if this condition is detected before 452c taking any steps, then an illegal input return 453c (istate = -3) occurs instead.) 454c -3 means illegal input was detected, before taking any 455c integration steps. see written message for details. 456c note.. if the solver detects an infinite loop of calls 457c to the solver with illegal input, it will cause 458c the run to stop. 459c -4 means there were repeated error test failures on 460c one attempted step, before completing the requested 461c task, but the integration was successful as far as t. 462c the problem may have a singularity, or the input 463c may be inappropriate. 464c -5 means there were repeated convergence test failures on 465c one attempted step, before completing the requested 466c task, but the integration was successful as far as t. 467c this may be caused by an inaccurate jacobian matrix, 468c if one is being used. 469c -6 means ewt(i) became zero for some i during the 470c integration. pure relative error control (atol(i)=0.0) 471c was requested on a variable which has now vanished. 472c the integration was successful as far as t. 473c 474c note.. since the normal output value of istate is 2, 475c it does not need to be reset for normal continuation. 476c also, since a negative input value of istate will be 477c regarded as illegal, a negative output value requires the 478c user to change it, and possibly other inputs, before 479c calling the solver again. 480c 481c iopt = an integer flag to specify whether or not any optional 482c inputs are being used on this call. input only. 483c the optional inputs are listed separately below. 484c iopt = 0 means no optional inputs are being used. 485c default values will be used in all cases. 486c iopt = 1 means one or more optional inputs are being used. 487c 488c rwork = a real working array (double precision). 489c the length of rwork must be at least 490c 20 + nyh*(maxord + 1) + 3*neq + lwm where 491c nyh = the initial value of neq, 492c maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a 493c smaller value is given as an optional input), 494c lwm = 0 if miter = 0, 495c lwm = neq**2 + 2 if miter is 1 or 2, 496c lwm = neq + 2 if miter = 3, and 497c lwm = (2*ml+mu+1)*neq + 2 if miter is 4 or 5. 498c (see the mf description for meth and miter.) 499c thus if maxord has its default value and neq is constant, 500c this length is.. 501c 20 + 16*neq for mf = 10, 502c 22 + 16*neq + neq**2 for mf = 11 or 12, 503c 22 + 17*neq for mf = 13, 504c 22 + 17*neq + (2*ml+mu)*neq for mf = 14 or 15, 505c 20 + 9*neq for mf = 20, 506c 22 + 9*neq + neq**2 for mf = 21 or 22, 507c 22 + 10*neq for mf = 23, 508c 22 + 10*neq + (2*ml+mu)*neq for mf = 24 or 25. 509c the first 20 words of rwork are reserved for conditional 510c and optional inputs and optional outputs. 511c 512c the following word in rwork is a conditional input.. 513c rwork(1) = tcrit = critical value of t which the solver 514c is not to overshoot. required if itask is 515c 4 or 5, and ignored otherwise. (see itask.) 516c 517c lrw = the length of the array rwork, as declared by the user. 518c (this will be checked by the solver.) 519c 520c iwork = an integer work array. the length of iwork must be at least 521c 20 if miter = 0 or 3 (mf = 10, 13, 20, 23), or 522c 20 + neq otherwise (mf = 11, 12, 14, 15, 21, 22, 24, 25). 523c the first few words of iwork are used for conditional and 524c optional inputs and optional outputs. 525c 526c the following 2 words in iwork are conditional inputs.. 527c iwork(1) = ml these are the lower and upper 528c iwork(2) = mu half-bandwidths, respectively, of the 529c banded jacobian, excluding the main diagonal. 530c the band is defined by the matrix locations 531c (i,j) with i-ml .le. j .le. i+mu. ml and mu 532c must satisfy 0 .le. ml,mu .le. neq-1. 533c these are required if miter is 4 or 5, and 534c ignored otherwise. ml and mu may in fact be 535c the band parameters for a matrix to which 536c df/dy is only approximately equal. 537c 538c liw = the length of the array iwork, as declared by the user. 539c (this will be checked by the solver.) 540c 541c note.. the work arrays must not be altered between calls to lsode 542c for the same problem, except possibly for the conditional and 543c optional inputs, and except for the last 3*neq words of rwork. 544c the latter space is used for internal scratch space, and so is 545c available for use by the user outside lsode between calls, if 546c desired (but not for use by f or jac). 547c 548c jac = the name of the user-supplied routine (miter = 1 or 4) to 549c compute the jacobian matrix, df/dy, as a function of 550c the scalar t and the vector y. it is to have the form 551c subroutine jac (neq, t, y, ml, mu, pd, nrowpd) 552c dimension y(*), pd(nrowpd,*) 553c where neq, t, y, ml, mu, and nrowpd are input and the array 554c pd is to be loaded with partial derivatives (elements of 555c the jacobian matrix) on output. pd must be given a first 556c dimension of nrowpd. t and y have the same meaning as in 557c subroutine f. (in the dimension statement above, 1 is a 558c dummy dimension.. it can be replaced by any value.) 559c in the full matrix case (miter = 1), ml and mu are 560c ignored, and the jacobian is to be loaded into pd in 561c columnwise manner, with df(i)/dy(j) loaded into pd(i,j). 562c in the band matrix case (miter = 4), the elements 563c within the band are to be loaded into pd in columnwise 564c manner, with diagonal lines of df/dy loaded into the rows 565c of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j). 566c ml and mu are the half-bandwidth parameters (see iwork). 567c the locations in pd in the two triangular areas which 568c correspond to nonexistent matrix elements can be ignored 569c or loaded arbitrarily, as they are overwritten by lsode. 570c jac need not provide df/dy exactly. a crude 571c approximation (possibly with a smaller bandwidth) will do. 572c in either case, pd is preset to zero by the solver, 573c so that only the nonzero elements need be loaded by jac. 574c each call to jac is preceded by a call to f with the same 575c arguments neq, t, and y. thus to gain some efficiency, 576c intermediate quantities shared by both calculations may be 577c saved in a user common block by f and not recomputed by jac, 578c if desired. also, jac may alter the y array, if desired. 579c jac must be declared external in the calling program. 580c subroutine jac may access user-defined quantities in 581c neq(2),... and y(neq(1)+1),... if neq is an array 582c (dimensioned in jac) and y has length exceeding neq(1). 583c see the descriptions of neq and y above. 584c 585c mf = the method flag. used only for input. the legal values of 586c mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, and 25. 587c mf has decimal digits meth and miter.. mf = 10*meth + miter. 588c meth indicates the basic linear multistep method.. 589c meth = 1 means the implicit adams method. 590c meth = 2 means the method based on backward 591c differentiation formulas (bdf-s). 592c miter indicates the corrector iteration method.. 593c miter = 0 means functional iteration (no jacobian matrix 594c is involved). 595c miter = 1 means chord iteration with a user-supplied 596c full (neq by neq) jacobian. 597c miter = 2 means chord iteration with an internally 598c generated (difference quotient) full jacobian 599c (using neq extra calls to f per df/dy value). 600c miter = 3 means chord iteration with an internally 601c generated diagonal jacobian approximation. 602c (using 1 extra call to f per df/dy evaluation). 603c miter = 4 means chord iteration with a user-supplied 604c banded jacobian. 605c miter = 5 means chord iteration with an internally 606c generated banded jacobian (using ml+mu+1 extra 607c calls to f per df/dy evaluation). 608c if miter = 1 or 4, the user must supply a subroutine jac 609c (the name is arbitrary) as described above under jac. 610c for other values of miter, a dummy argument can be used. 611c 612c!optional inputs. 613c 614c the following is a list of the optional inputs provided for in the 615c call sequence. (see also part ii.) for each such input variable, 616c this table lists its name as used in this documentation, its 617c location in the call sequence, its meaning, and the default value. 618c the use of any of these inputs requires iopt = 1, and in that 619c case all of these inputs are examined. a value of zero for any 620c of these optional inputs will cause the default value to be used. 621c thus to use a subset of the optional inputs, simply preload 622c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and 623c then set those of interest to nonzero values. 624c 625c name location meaning and default value 626c 627c h0 rwork(5) the step size to be attempted on the first step. 628c the default value is determined by the solver. 629c 630c hmax rwork(6) the maximum absolute step size allowed. 631c the default value is infinite. 632c 633c hmin rwork(7) the minimum absolute step size allowed. 634c the default value is 0. (this lower bound is not 635c enforced on the final step before reaching tcrit 636c when itask = 4 or 5.) 637c 638c maxord iwork(5) the maximum order to be allowed. the default 639c value is 12 if meth = 1, and 5 if meth = 2. 640c if maxord exceeds the default value, it will 641c be reduced to the default value. 642c if maxord is changed during the problem, it may 643c cause the current order to be reduced. 644c 645c mxstep iwork(6) maximum number of (internally defined) steps 646c allowed during one call to the solver. 647c the default value is 500. 648c 649c mxhnil iwork(7) maximum number of messages printed (per problem) 650c warning that t + h = t on a step (h = step size). 651c this must be positive to result in a non-default 652c value. the default value is 10. 653c 654c!optional outputs. 655c 656c as optional additional output from lsode, the variables listed 657c below are quantities related to the performance of lsode 658c which are available to the user. these are communicated by way of 659c the work arrays, but also have internal mnemonic names as shown. 660c except where stated otherwise, all of these outputs are defined 661c on any successful return from lsode, and on any return with 662c istate = -1, -2, -4, -5, or -6. on an illegal input return 663c (istate = -3), they will be unchanged from their existing values 664c (if any), except possibly for tolsf, lenrw, and leniw. 665c on any error return, outputs relevant to the error will be defined, 666c as noted below. 667c 668c name location meaning 669c 670c hu rwork(11) the step size in t last used (successfully). 671c 672c hcur rwork(12) the step size to be attempted on the next step. 673c 674c tcur rwork(13) the current value of the independent variable 675c which the solver has actually reached, i.e. the 676c current internal mesh point in t. on output, tcur 677c will always be at least as far as the argument 678c t, but may be farther (if interpolation was done). 679c 680c tolsf rwork(14) a tolerance scale factor, greater than 1.0, 681c computed when a request for too much accuracy was 682c detected (istate = -3 if detected at the start of 683c the problem, istate = -2 otherwise). if itol is 684c left unaltered but rtol and atol are uniformly 685c scaled up by a factor of tolsf for the next call, 686c then the solver is deemed likely to succeed. 687c (the user may also ignore tolsf and alter the 688c tolerance parameters in any other way appropriate.) 689c 690c nst iwork(11) the number of steps taken for the problem so far. 691c 692c nfe iwork(12) the number of f evaluations for the problem so far. 693c 694c nje iwork(13) the number of jacobian evaluations (and of matrix 695c lu decompositions) for the problem so far. 696c 697c nqu iwork(14) the method order last used (successfully). 698c 699c nqcur iwork(15) the order to be attempted on the next step. 700c 701c imxer iwork(16) the index of the component of largest magnitude in 702c the weighted local error vector ( e(i)/ewt(i) ), 703c on an error return with istate = -4 or -5. 704c 705c lenrw iwork(17) the length of rwork actually required. 706c this is defined on normal returns and on an illegal 707c input return for insufficient storage. 708c 709c leniw iwork(18) the length of iwork actually required. 710c this is defined on normal returns and on an illegal 711c input return for insufficient storage. 712c 713c the following two arrays are segments of the rwork array which 714c may also be of interest to the user as optional outputs. 715c for each array, the table below gives its internal name, 716c its base address in rwork, and its description. 717c 718c name base address description 719c 720c yh 21 the nordsieck history array, of size nyh by 721c (nqcur + 1), where nyh is the initial value 722c of neq. for j = 0,1,...,nqcur, column j+1 723c of yh contains hcur**j/factorial(j) times 724c the j-th derivative of the interpolating 725c polynomial currently representing the solution, 726c evaluated at t = tcur. 727c 728c acor lenrw-neq+1 array of size neq used for the accumulated 729c corrections on each step, scaled on output 730c to represent the estimated local error in y 731c on the last step. this is the vector e in 732c the description of the error control. it is 733c defined only on a successful return from lsode. 734c 735c 736c!part ii. other routines callable. 737c 738c the following are optional calls which the user may make to 739c gain additional capabilities in conjunction with lsode. 740c (the routines xsetun and xsetf are designed to conform to the 741c slatec error handling package.) 742c 743c form of call function 744c call xsetun(lun) set the logical unit number, lun, for 745c output of messages from lsode, if 746c the default is not desired. 747c the default value of lun is 6. 748c 749c call xsetf(mflag) set a flag to control the printing of 750c messages by lsode. 751c mflag = 0 means do not print. (danger.. 752c this risks losing valuable information.) 753c mflag = 1 means print (the default). 754c 755c either of the above calls may be made at 756c any time and will take effect immediately. 757c 758c call svcom (rsav, isav) store in rsav and isav the contents 759c of the internal common blocks used by 760c lsode (see part iii below). 761c rsav must be a real array of length 219 762c or more, and isav must be an integer 763c array of length 41 or more. 764c 765c call rscom (rsav, isav) restore, from rsav and isav, the contents 766c of the internal common blocks used by 767c lsode. presumes a prior call to svcom 768c with the same arguments. 769c 770c svcom and rscom are useful if 771c interrupting a run and restarting 772c later, or alternating between two or 773c more problems solved with lsode. 774c 775c call intdy(,,,,,) provide derivatives of y, of various 776c (see below) orders, at a specified point t, if 777c desired. it may be called only after 778c a successful return from lsode. 779c 780c the detailed instructions for using intdy are as follows. 781c the form of the call is.. 782c 783c call intdy (t, k, rwork(21), nyh, dky, iflag) 784c 785c the input parameters are.. 786c 787c t = value of independent variable where answers are desired 788c (normally the same as the t last returned by lsode). 789c for valid results, t must lie between tcur - hu and tcur. 790c (see optional outputs for tcur and hu.) 791c k = integer order of the derivative desired. k must satisfy 792c 0 .le. k .le. nqcur, where nqcur is the current order 793c (see optional outputs). the capability corresponding 794c to k = 0, i.e. computing y(t), is already provided 795c by lsode directly. since nqcur .ge. 1, the first 796c derivative dy/dt is always available with intdy. 797c rwork(21) = the base address of the history array yh. 798c nyh = column length of yh, equal to the initial value of neq. 799c 800c the output parameters are.. 801c 802c dky = a real array of length neq containing the computed value 803c of the k-th derivative of y(t). 804c iflag = integer flag, returned as 0 if k and t were legal, 805c -1 if k was illegal, and -2 if t was illegal. 806c on an error return, a message is also written. 807c 808c!part iii. common blocks. 809c 810c if lsode is to be used in an overlay situation, the user 811c must declare, in the primary overlay, the variables in.. 812c (1) the call sequence to lsode, 813c (2) the two internal common blocks 814c /ls0001/ of length 258 (219 double precision words 815c followed by 39 integer words), 816c /eh0001/ of length 2 (integer words). 817c 818c if lsode is used on a system in which the contents of internal 819c common blocks are not preserved between calls, the user should 820c declare the above two common blocks in his main program to insure 821c that their contents are preserved. 822c 823c if the solution of a given problem by lsode is to be interrupted 824c and then later continued, such as when restarting an interrupted run 825c or alternating between two or more problems, the user should save, 826c following the return from the last lsode call prior to the 827c interruption, the contents of the call sequence variables and the 828c internal common blocks, and later restore these values before the 829c next lsode call for that problem. to save and restore the common 830c blocks, use subroutines svcom and rscom (see part ii above). 831c 832c 833c!part iv. optionally replaceable solver routines. 834c 835c below are descriptions of two routines in the lsode package which 836c relate to the measurement of errors. either routine can be 837c replaced by a user-supplied version, if desired. however, since such 838c a replacement may have a major impact on performance, it should be 839c done only when absolutely necessary, and only with great caution. 840c (note.. the means by which the package version of a routine is 841c superseded by the user-s version may be system-dependent.) 842c 843c (a) ewset. 844c the following subroutine is called just before each internal 845c integration step, and sets the array of error weights, ewt, as 846c described under itol/rtol/atol above.. 847c subroutine ewset (neq, itol, rtol, atol, ycur, ewt) 848c where neq, itol, rtol, and atol are as in the lsode call sequence, 849c ycur contains the current dependent variable vector, and 850c ewt is the array of weights set by ewset. 851c 852c if the user supplies this subroutine, it must return in ewt(i) 853c (i = 1,...,neq) a positive quantity suitable for comparing errors 854c in y(i) to. the ewt array returned by ewset is passed to the 855c vnorm routine (see below), and also used by lsode in the computation 856c of the optional output imxer, the diagonal jacobian approximation, 857c and the increments for difference quotient jacobians. 858c 859c in the user-supplied version of ewset, it may be desirable to use 860c the current values of derivatives of y. derivatives up to order nq 861c are available from the history array yh, described above under 862c optional outputs. in ewset, yh is identical to the ycur array, 863c extended to nq + 1 columns with a column length of nyh and scale 864c factors of h**j/factorial(j). on the first call for the problem, 865c given by nst = 0, nq is 1 and h is temporarily set to 1.0. 866c the quantities nq, nyh, h, and nst can be obtained by including 867c in ewset the statements.. 868c double precision h, rls 869c common /ls0001/ rls(219),ils(39) 870c nq = ils(35) 871c nyh = ils(14) 872c nst = ils(36) 873c h = rls(213) 874c thus, for example, the current value of dy/dt can be obtained as 875c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is 876c unnecessary when nst = 0). 877c 878c (b) vnorm. 879c the following is a real function routine which computes the weighted 880c root-mean-square norm of a vector v.. 881c d = vnorm (n, v, w) 882c where.. 883c n = the length of the vector, 884c v = real array of length n containing the vector, 885c w = real array of length n containing weights, 886c d = sqrt( (1/n) * sum(v(i)*w(i))**2 ). 887c vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where 888c ewt is as set by subroutine ewset. 889c 890c if the user supplies this function, it should return a non-negative 891c value of vnorm suitable for use in the error control in lsode. 892c none of the arguments should be altered by vnorm. 893c for example, a user-supplied vnorm routine might.. 894c -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or 895c -ignore some components of v in the norm, with the effect of 896c suppressing the error control on those components of y. 897c 898c 899c!other routines in the lsode package. 900c 901c in addition to subroutine lsode, the lsode package includes the 902c following subroutines and function routines.. 903c intdy computes an interpolated value of the y vector at t = tout. 904c stode is the core integrator, which does one step of the 905c integration and the associated error control. 906c cfode sets all method coefficients and test constants. 907c prepj computes and preprocesses the jacobian matrix j = df/dy 908c and the newton iteration matrix p = i - h*l0*j. 909c solsy manages solution of linear system in chord iteration. 910c ewset sets the error weight vector ewt before each step. 911c vnorm computes the weighted r.m.s. norm of a vector. 912c svcom and rscom are user-callable routines to save and restore, 913c respectively, the contents of the internal common blocks. 914c dgefa and dgesl are routines from linpack for solving full 915c systems of linear algebraic equations. 916c dgbfa and dgbsl are routines from linpack for solving banded 917c linear systems. 918c daxpy, dscal, idamax, and ddot are basic linear algebra modules 919c (blas) used by the above linpack routines. 920c dlamch computes the unit roundoff in a machine-independent manner. 921c xerrwv, xsetun, and xsetf handle the printing of all error 922c messages and warnings. xerrwv is machine-dependent. 923c note.. vnorm, idamax, ddot, and dlamch are function routines. 924c all the others are subroutines. 925c 926c the intrinsic and external routines used by lsode are.. 927c abs, max, min, dble, max, min, mod, sign, sqrt, and write. 928c 929c a block data subprogram is also included with the package, 930c for loading some of the variables in internal common. 931c 932c!author and contact 933c alan c. hindmarsh, 934c mathematics and statistics division, l-316 935c lawrence livermore national laboratory 936c livermore, ca 94550. 937c 938c!reference. 939c alan c. hindmarsh, lsode and lsodi, two new initial value 940c ordinary differential equation solvers, 941c acm-signum newsletter, vol. 15, no. 4 (1980), pp. 10-11. 942c 943c! 944c this is the august 13, 1981 version of lsode. 945c----------------------------------------------------------------------- 946c the following card is for optimized compilation on llnl compilers. 947clll. optimize 948c----------------------------------------------------------------------- 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 tret, 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 dlamch, vnorm 961 dimension mord(2) 962 logical ihit 963c----------------------------------------------------------------------- 964c the following internal common block contains 965c (a) variables which are local to any subroutine but whose values must 966c be preserved between calls to the routine (own variables), and 967c (b) variables which are communicated between subroutines. 968c the structure of the block is as follows.. all real variables are 969c listed first, followed by all integers. within each type, the 970c variables are grouped with those local to subroutine lsode first, 971c then those local to subroutine stode, and finally those used 972c for communication. the block is declared in subroutines 973c lsode, intdy, stode, prepj, and solsy. groups of variables are 974c replaced by dummy arrays in the common declarations in routines 975c where those variables are not used. 976c----------------------------------------------------------------------- 977cDEC$ ATTRIBUTES DLLIMPORT:: /ls0001/ 978 common /ls0001/ tret, rowns(209), 979 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 980 2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, 981 3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), 982 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 983 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu 984 save/ls0001/ 985c 986 data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/ 987c----------------------------------------------------------------------- 988c block a. 989c this code block is executed on every call. 990c it tests istate and itask for legality and branches appropiately. 991c if istate .gt. 1 but the flag init shows that initialization has 992c not yet been done, an error return occurs. 993c if istate = 1 and tout = t, jump to block g and return immediately. 994c----------------------------------------------------------------------- 995 ierror=0 996 if (istate .lt. 1 .or. istate .gt. 3) go to 601 997 if (itask .lt. 1 .or. itask .gt. 5) go to 602 998 if (istate .eq. 1) go to 10 999 if (init .eq. 0) go to 603 1000 if (istate .eq. 2) go to 200 1001 go to 20 1002 10 init = 0 1003 if (tout .eq. t) go to 430 1004 20 ntrep = 0 1005c----------------------------------------------------------------------- 1006c block b. 1007c the next code block is executed for the initial call (istate = 1), 1008c or for a continuation call with parameter changes (istate = 3). 1009c it contains checking of all inputs and various initializations. 1010c 1011c first check legality of the non-optional inputs neq, itol, iopt, 1012c mf, ml, and mu. 1013c----------------------------------------------------------------------- 1014 if (neq(1) .le. 0) go to 604 1015 if (istate .eq. 1) go to 25 1016 if (neq(1) .gt. n) go to 605 1017 25 n = neq(1) 1018 if (itol .lt. 1 .or. itol .gt. 4) go to 606 1019 if (iopt .lt. 0 .or. iopt .gt. 1) go to 607 1020 meth = mf/10 1021 miter = mf - 10*meth 1022 if (meth .lt. 1 .or. meth .gt. 2) go to 608 1023 if (miter .lt. 0 .or. miter .gt. 5) go to 608 1024 if (miter .le. 3) go to 30 1025 ml = iwork(1) 1026 mu = iwork(2) 1027 if (ml .lt. 0 .or. ml .ge. n) go to 609 1028 if (mu .lt. 0 .or. mu .ge. n) go to 610 1029 30 continue 1030c next process and check the optional inputs. -------------------------- 1031 if (iopt .eq. 1) go to 40 1032 maxord = mord(meth) 1033 mxstep = mxstp0 1034 mxhnil = mxhnl0 1035 if (istate .eq. 1) h0 = 0.0d+0 1036 hmxi = 0.0d+0 1037 hmin = 0.0d+0 1038 go to 60 1039 40 maxord = iwork(5) 1040 if (maxord .lt. 0) go to 611 1041 if (maxord .eq. 0) maxord = 100 1042 maxord = min(maxord,mord(meth)) 1043 mxstep = iwork(6) 1044 if (mxstep .lt. 0) go to 612 1045 if (mxstep .eq. 0) mxstep = mxstp0 1046 mxhnil = iwork(7) 1047 if (mxhnil .lt. 0) go to 613 1048 if (mxhnil .eq. 0) mxhnil = mxhnl0 1049 if (istate .ne. 1) go to 50 1050 h0 = rwork(5) 1051 if ((tout - t)*h0 .lt. 0.0d+0) go to 614 1052 50 hmax = rwork(6) 1053 if (hmax .lt. 0.0d+0) go to 615 1054 hmxi = 0.0d+0 1055 if (hmax .gt. 0.0d+0) hmxi = 1.0d+0/hmax 1056 hmin = rwork(7) 1057 if (hmin .lt. 0.0d+0) go to 616 1058c----------------------------------------------------------------------- 1059c set work array pointers and check lengths lrw and liw. 1060c pointers to segments of rwork and iwork are named by prefixing l to 1061c the name of the segment. e.g., the segment yh starts at rwork(lyh). 1062c segments of rwork (in order) are denoted yh, wm, ewt, savf, acor. 1063c----------------------------------------------------------------------- 1064 60 lyh = 21 1065 if (istate .eq. 1) nyh = n 1066 lwm = lyh + (maxord + 1)*nyh 1067 if (miter .eq. 0) lenwm = 0 1068 if (miter .eq. 1 .or. miter .eq. 2) lenwm = n*n + 2 1069 if (miter .eq. 3) lenwm = n + 2 1070 if (miter .ge. 4) lenwm = (2*ml + mu + 1)*n + 2 1071 lewt = lwm + lenwm 1072 lsavf = lewt + n 1073 lacor = lsavf + n 1074 lenrw = lacor + n - 1 1075 iwork(17) = lenrw 1076 liwm = 1 1077 leniw = 20 + n 1078 if (miter .eq. 0 .or. miter .eq. 3) leniw = 20 1079 iwork(18) = leniw 1080 if (lenrw .gt. lrw) go to 617 1081 if (leniw .gt. liw) go to 618 1082c check rtol and atol for legality. ------------------------------------ 1083 rtoli = rtol(1) 1084 atoli = atol(1) 1085 do 70 i = 1,n 1086 if (itol .ge. 3) rtoli = rtol(i) 1087 if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) 1088 if (rtoli .lt. 0.0d+0) go to 619 1089 if (atoli .lt. 0.0d+0) go to 620 1090 70 continue 1091 if (istate .eq. 1) go to 100 1092c if istate = 3, set flag to signal parameter changes to stode. -------- 1093 jstart = -1 1094 if (nq .le. maxord) go to 90 1095c maxord was reduced below nq. copy yh(*,maxord+2) into savf. --------- 1096 do 80 i = 1,n 1097 80 rwork(i+lsavf-1) = rwork(i+lwm-1) 1098c reload wm(1) = rwork(lwm), since lwm may have changed. --------------- 1099 90 if (miter .gt. 0) rwork(lwm) = sqrt(uround) 1100 if (n .eq. nyh) go to 200 1101c neq was reduced. zero part of yh to avoid undefined references. ----- 1102 i1 = lyh + l*nyh 1103 i2 = lyh + (maxord + 1)*nyh - 1 1104 if (i1 .gt. i2) go to 200 1105 do 95 i = i1,i2 1106 95 rwork(i) = 0.0d+0 1107 go to 200 1108c----------------------------------------------------------------------- 1109c block c. 1110c the next block is for the initial call only (istate = 1). 1111c it contains all remaining initializations, the initial call to f, 1112c and the calculation of the initial step size. 1113c the error weights in ewt are inverted after being loaded. 1114c----------------------------------------------------------------------- 1115 100 uround = dlamch('p') 1116 tn = t 1117 if (itask .ne. 4 .and. itask .ne. 5) go to 110 1118 tcrit = rwork(1) 1119 if ((tcrit - tout)*(tout - t) .lt. 0.0d+0) go to 625 1120 if (h0 .ne. 0.0d+0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d+0) 1121 1 h0 = tcrit - t 1122 110 jstart = 0 1123 if (miter .gt. 0) rwork(lwm) = sqrt(uround) 1124 nhnil = 0 1125 nst = 0 1126 nje = 0 1127 nslast = 0 1128 hu = 0.0d+0 1129 nqu = 0 1130 ccmax = 0.30d+0 1131 maxcor = 3 1132 msbp = 20 1133 mxncf = 10 1134c initial call to f. (lf0 points to yh(*,2).) ------------------------- 1135 lf0 = lyh + nyh 1136 call f (neq, t, y, rwork(lf0)) 1137 if(ierror.gt.0) return 1138 nfe = 1 1139c load the initial value vector in yh. --------------------------------- 1140 do 115 i = 1,n 1141 115 rwork(i+lyh-1) = y(i) 1142c load and invert the ewt array. (h is temporarily set to 1.0.) ------- 1143 nq = 1 1144 h = 1.0d+0 1145 call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) 1146 do 120 i = 1,n 1147 if (rwork(i+lewt-1) .le. 0.0d+0) go to 621 1148 120 rwork(i+lewt-1) = 1.0d+0/rwork(i+lewt-1) 1149c----------------------------------------------------------------------- 1150c the coding below computes the step size, h0, to be attempted on the 1151c first step, unless the user has supplied a value for this. 1152c first check that tout - t differs significantly from zero. 1153c a scalar tolerance quantity tol is computed, as max(rtol(i)) 1154c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted 1155c so as to be between 100*uround and 1.0e-3. 1156c then the computed value h0 is given by.. 1157c neq 1158c h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2 ) 1159c 1 1160c where w0 = max ( abs(t), abs(tout) ), 1161c f(i) = i-th component of initial value of f, 1162c ywt(i) = ewt(i)/tol (a weight for y(i)). 1163c the sign of h0 is inferred from the initial values of tout and t. 1164c----------------------------------------------------------------------- 1165 if (h0 .ne. 0.0d+0) go to 180 1166 tdist = abs(tout - t) 1167 w0 = max(abs(t),abs(tout)) 1168 if (tdist .lt. 2.0d+0*uround*w0) go to 622 1169 tol = rtol(1) 1170 if (itol .le. 2) go to 140 1171 do 130 i = 1,n 1172 130 tol = max(tol,rtol(i)) 1173 140 if (tol .gt. 0.0d+0) go to 160 1174 atoli = atol(1) 1175 do 150 i = 1,n 1176 if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i) 1177 ayi = abs(y(i)) 1178 if (ayi .ne. 0.0d+0) tol = max(tol,atoli/ayi) 1179 150 continue 1180 160 tol = max(tol,100.0d+0*uround) 1181 tol = min(tol,0.0010d+0) 1182 sum = vnorm (n, rwork(lf0), rwork(lewt)) 1183 sum = 1.0d+0/(tol*w0*w0) + tol*sum**2 1184 h0 = 1.0d+0/sqrt(sum) 1185 h0 = min(h0,tdist) 1186 h0 = sign(h0,tout-t) 1187c adjust h0 if necessary to meet hmax bound. --------------------------- 1188 180 rh = abs(h0)*hmxi 1189 if (rh .gt. 1.0d+0) h0 = h0/rh 1190c load h with h0 and scale yh(*,2) by h0. ------------------------------ 1191 h = h0 1192 do 190 i = 1,n 1193 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1) 1194 go to 270 1195c----------------------------------------------------------------------- 1196c block d. 1197c the next code block is for continuation calls only (istate = 2 or 3) 1198c and is to check stop conditions before taking a step. 1199c----------------------------------------------------------------------- 1200 200 nslast = nst 1201 go to (210, 250, 220, 230, 240), itask 1202 210 if ((tn - tout)*h .lt. 0.0d+0) go to 250 1203 call intdy (tout, 0, rwork(lyh), nyh, y, iflag) 1204 if (iflag .ne. 0) go to 627 1205 t = tout 1206 go to 420 1207 220 tp = tn - hu*(1.0d+0 + 100.0d+0*uround) 1208 if ((tp - tout)*h .gt. 0.0d+0) go to 623 1209 if ((tn - tout)*h .lt. 0.0d+0) go to 250 1210 go to 400 1211 230 tcrit = rwork(1) 1212 if ((tn - tcrit)*h .gt. 0.0d+0) go to 624 1213 if ((tcrit - tout)*h .lt. 0.0d+0) go to 625 1214 if ((tn - tout)*h .lt. 0.0d+0) go to 245 1215 call intdy (tout, 0, rwork(lyh), nyh, y, iflag) 1216 if (iflag .ne. 0) go to 627 1217 t = tout 1218 go to 420 1219 240 tcrit = rwork(1) 1220 if ((tn - tcrit)*h .gt. 0.0d+0) go to 624 1221 245 hmx = abs(tn) + abs(h) 1222 ihit = abs(tn - tcrit) .le. 100.0d+0*uround*hmx 1223 if (ihit) go to 400 1224 tnext = tn + h*(1.0d+0 + 4.0d+0*uround) 1225 if ((tnext - tcrit)*h .le. 0.0d+0) go to 250 1226 h = (tcrit - tn)*(1.0d+0 - 4.0d+0*uround) 1227 if (istate .eq. 2) jstart = -2 1228c----------------------------------------------------------------------- 1229c block e. 1230c the next block is normally executed for all calls and contains 1231c the call to the one-step core integrator stode. 1232c 1233c this is a looping point for the integration steps. 1234c 1235c first check for too many steps being taken, update ewt (if not at 1236c start of problem), check for too much accuracy being requested, and 1237c check for h below the roundoff level in t. 1238c----------------------------------------------------------------------- 1239 250 continue 1240 if ((nst-nslast) .ge. mxstep) go to 500 1241 call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt)) 1242 do 260 i = 1,n 1243 if (rwork(i+lewt-1) .le. 0.0d+0) go to 510 1244 260 rwork(i+lewt-1) = 1.0d+0/rwork(i+lewt-1) 1245 270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt)) 1246 if (tolsf .le. 1.0d+0) go to 280 1247 tolsf = tolsf*2.0d+0 1248 if (nst .eq. 0) go to 626 1249 go to 520 1250 280 if ((tn + h) .ne. tn) go to 290 1251 nhnil = nhnil + 1 1252 if (nhnil .gt. mxhnil) go to 290 1253 call xerrwv('lsode-- caution... t (=r1) and h (=r2) are', 1254 1 50, 101, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1255 call xerrwv( 1256 1 ' such that t + h = t at next step', 1257 1 60, 101, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1258 call xerrwv(' (h = step). integration continues', 1259 1 50, 101, 1, 0, 0, 0, 2, tn, h) 1260 if (nhnil .lt. mxhnil) go to 290 1261 call xerrwv('lsode-- preceding message given i1 times', 1262 1 50, 102, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1263 call xerrwv(' will not be repeated', 1264 1 50, 102, 1, 1, mxhnil, 0, 0, 0.0d+0, 0.0d+0) 1265 290 continue 1266c----------------------------------------------------------------------- 1267c call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,jac,prepj,solsy) 1268c----------------------------------------------------------------------- 1269 call stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt), 1270 1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm), 1271 2 f, jac, prepj, solsy) 1272 if(ierror.gt.0) return 1273 kgo = 1 - kflag 1274 go to (300, 530, 540), kgo 1275c----------------------------------------------------------------------- 1276c block f. 1277c the following block handles the case of a successful return from the 1278c core integrator (kflag = 0). test for stop conditions. 1279c----------------------------------------------------------------------- 1280 300 init = 1 1281 go to (310, 400, 330, 340, 350), itask 1282c itask = 1. if tout has been reached, interpolate. ------------------- 1283 310 if ((tn - tout)*h .lt. 0.0d+0) go to 250 1284 call intdy (tout, 0, rwork(lyh), nyh, y, iflag) 1285 t = tout 1286 go to 420 1287c itask = 3. jump to exit if tout was reached. ------------------------ 1288 330 if ((tn - tout)*h .ge. 0.0d+0) go to 400 1289 go to 250 1290c itask = 4. see if tout or tcrit was reached. adjust h if necessary. 1291 340 if ((tn - tout)*h .lt. 0.0d+0) go to 345 1292 call intdy (tout, 0, rwork(lyh), nyh, y, iflag) 1293 t = tout 1294 go to 420 1295 345 hmx = abs(tn) + abs(h) 1296 ihit = abs(tn - tcrit) .le. 100.0d+0*uround*hmx 1297 if (ihit) go to 400 1298 tnext = tn + h*(1.0d+0 + 4.0d+0*uround) 1299 if ((tnext - tcrit)*h .le. 0.0d+0) go to 250 1300 h = (tcrit - tn)*(1.0d+0 - 4.0d+0*uround) 1301 jstart = -2 1302 go to 250 1303c itask = 5. see if tcrit was reached and jump to exit. --------------- 1304 350 hmx = abs(tn) + abs(h) 1305 ihit = abs(tn - tcrit) .le. 100.0d+0*uround*hmx 1306c----------------------------------------------------------------------- 1307c block g. 1308c the following block handles all successful returns from lsode. 1309c if itask .ne. 1, y is loaded from yh and t is set accordingly. 1310c istate is set to 2, the illegal input counter is zeroed, and the 1311c optional outputs are loaded into the work arrays before returning. 1312c if istate = 1 and tout = t, there is a return with no action taken, 1313c except that if this has happened repeatedly, the run is terminated. 1314c----------------------------------------------------------------------- 1315 400 do 410 i = 1,n 1316 410 y(i) = rwork(i+lyh-1) 1317 t = tn 1318 if (itask .ne. 4 .and. itask .ne. 5) go to 420 1319 if (ihit) t = tcrit 1320 420 istate = 2 1321 illin = 0 1322 rwork(11) = hu 1323 rwork(12) = h 1324 rwork(13) = tn 1325 iwork(11) = nst 1326 iwork(12) = nfe 1327 iwork(13) = nje 1328 iwork(14) = nqu 1329 iwork(15) = nq 1330 return 1331c 1332 430 ntrep = ntrep + 1 1333 if (ntrep .lt. 5) return 1334 call xerrwv( 1335 1 'lsode-- calls with istate = 1 and tout = t (=r1) ', 1336 1 60, 301, 1, 0, 0, 0, 1, t, 0.0d+0) 1337 go to 800 1338c----------------------------------------------------------------------- 1339c block h. 1340c the following block handles all unsuccessful returns other than 1341c those for illegal input. first the error message routine is called. 1342c if there was an error test or convergence test failure, imxer is set. 1343c then y is loaded from yh, t is set to tn, and the illegal input 1344c counter illin is set to 0. the optional outputs are loaded into 1345c the work arrays before returning. 1346c----------------------------------------------------------------------- 1347c the maximum number of steps was taken before reaching tout. ---------- 1348 500 call xerrwv('lsode-- at t (=r1), mxstep (=i1) steps ', 1349 1 50, 201, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1350 call xerrwv('necessary before reaching tout', 1351 1 50, 201, 1, 1, mxstep, 0, 1, tn, 0.0d+0) 1352 istate = -1 1353 go to 580 1354c ewt(i) .le. 0.0 for some i (not at start of problem). ---------------- 1355 510 ewti = rwork(lewt+i-1) 1356 call xerrwv('lsode-- at t (=r1),ewt(i1) (=r2) is .le.0', 1357 1 50, 202, 1, 1, i, 0, 2, tn, ewti) 1358 istate = -6 1359 go to 580 1360c too much accuracy requested for machine precision. ------------------- 1361 520 call xerrwv('lsode-- a t (=r1), too much precision required', 1362 1 50, 203, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1363 call xerrwv(' w.r.t. machine precision tolsf (=r2) ', 1364 1 50, 203, 1, 0, 0, 0, 2, tn, tolsf) 1365 rwork(14) = tolsf 1366 istate = -2 1367 go to 580 1368c kflag = -1. error test failed repeatedly or with abs(h) = hmin. ----- 1369 530 call xerrwv('lsode-- at t(=r1) for step h(=r2), error test', 1370 1 50, 204, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1371 call xerrwv(' failed with abs(h) = hmin', 1372 1 50, 204, 1, 0, 0, 0, 2, tn, h) 1373 istate = -4 1374 go to 560 1375c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ---- 1376 540 call xerrwv('lsode-- at t (=r1) with step h (=r2), ' , 1377 1 50, 205, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1378 call xerrwv(' corrector does not converge ', 1379 1 50, 205, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1380 call xerrwv(' with abs(h) = hmin ', 1381 1 30, 205, 1, 0, 0, 0, 2, tn, h) 1382 istate = -5 1383c compute imxer if relevant. ------------------------------------------- 1384 560 big = 0.0d+0 1385 imxer = 1 1386 do 570 i = 1,n 1387 size = abs(rwork(i+lacor-1)*rwork(i+lewt-1)) 1388 if (big .ge. size) go to 570 1389 big = size 1390 imxer = i 1391 570 continue 1392 iwork(16) = imxer 1393c set y vector, t, illin, and optional outputs. ------------------------ 1394 580 do 590 i = 1,n 1395 590 y(i) = rwork(i+lyh-1) 1396 t = tn 1397 illin = 0 1398 rwork(11) = hu 1399 rwork(12) = h 1400 rwork(13) = tn 1401 iwork(11) = nst 1402 iwork(12) = nfe 1403 iwork(13) = nje 1404 iwork(14) = nqu 1405 iwork(15) = nq 1406 return 1407c----------------------------------------------------------------------- 1408c block i. 1409c the following block handles all error returns due to illegal input 1410c (istate = -3), as detected before calling the core integrator. 1411c first the error message routine is called. then if there have been 1412c 5 consecutive such returns just before this call to the solver, 1413c the run is halted. 1414c----------------------------------------------------------------------- 1415 601 call xerrwv('lsode-- istate (=i1) illegal ', 1416 1 30, 1, 1, 1, istate, 0, 0, 0.0d+0, 0.0d+0) 1417 go to 700 1418 602 call xerrwv('lsode-- itask (=i1) illegal ', 1419 1 30, 2, 1, 1, itask, 0, 0, 0.0d+0, 0.0d+0) 1420 go to 700 1421 603 call xerrwv('lsode-- istate .gt. 1 ', 1422 1 50, 3, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1423 go to 700 1424 604 call xerrwv('lsode-- neq (=i1) .lt. 1 ', 1425 1 30, 4, 1, 1, neq(1), 0, 0, 0.0d+0, 0.0d+0) 1426 go to 700 1427 605 call xerrwv('lsode-- istate and neq increased from i1 to i2' , 1428 1 50, 5, 1, 2, n, neq(1), 0, 0.0d+0, 0.0d+0) 1429 go to 700 1430 606 call xerrwv('lsode-- itol (=i1) illegal ', 1431 1 30, 6, 1, 1, itol, 0, 0, 0.0d+0, 0.0d+0) 1432 go to 700 1433 607 call xerrwv('lsode-- iopt (=i1) illegal ', 1434 1 30, 7, 1, 1, iopt, 0, 0, 0.0d+0, 0.0d+0) 1435 go to 700 1436 608 call xerrwv('lsode-- mf (=i1) illegal ', 1437 1 30, 8, 1, 1, mf, 0, 0, 0.0d+0, 0.0d+0) 1438 go to 700 1439 609 call xerrwv('lsode-- ml (=i1) illegal.. .lt.0 or .ge.neq (=i2)', 1440 1 50, 9, 1, 2, ml, neq(1), 0, 0.0d+0, 0.0d+0) 1441 go to 700 1442 610 call xerrwv('lsode-- mu (=i1) illegal.. .lt.0 or .ge.neq (=i2)', 1443 1 50, 10, 1, 2, mu, neq(1), 0, 0.0d+0, 0.0d+0) 1444 go to 700 1445 611 call xerrwv('lsode-- maxord (=i1) .lt. 0 ', 1446 1 30, 11, 1, 1, maxord, 0, 0, 0.0d+0, 0.0d+0) 1447 go to 700 1448 612 call xerrwv('lsode-- mxstep (=i1) .lt. 0 ', 1449 1 30, 12, 1, 1, mxstep, 0, 0, 0.0d+0, 0.0d+0) 1450 go to 700 1451 613 call xerrwv('lsode-- mxhnil (=i1) .lt. 0 ', 1452 1 30, 13, 1, 1, mxhnil, 0, 0, 0.0d+0, 0.0d+0) 1453 go to 700 1454 614 call xerrwv('lsode-- tout (=r1) .gt. t (=r2) ', 1455 1 40, 14, 1, 0, 0, 0, 2, tout, t) 1456 call xerrwv(' h0 (=r1) gives integration direction' , 1457 1 50, 14, 1, 0, 0, 0, 1, h0, 0.0d+0) 1458 go to 700 1459 615 call xerrwv('lsode-- hmax (=r1) .lt. 0.0 ', 1460 1 30, 15, 1, 0, 0, 0, 1, hmax, 0.0d+0) 1461 go to 700 1462 616 call xerrwv('lsode-- hmin (=r1) .lt. 0.0 ', 1463 1 30, 16, 1, 0, 0, 0, 1, hmin, 0.0d+0) 1464 go to 700 1465 617 call xerrwv( 1466 1 'lsode-- necessary size for rwork (i1) larger than i2', 1467 1 60, 17, 1, 2, lenrw, lrw, 0, 0.0d+0, 0.0d+0) 1468 go to 700 1469 618 call xerrwv( 1470 1 'lsode-- necessary size for iwork (i1) larger than liw (i2)', 1471 1 60, 18, 1, 2, leniw, liw, 0, 0.0d+0, 0.0d+0) 1472 go to 700 1473 619 call xerrwv('lsode-- rtol(i1) is r1 .lt. 0.0 ', 1474 1 40, 19, 1, 1, i, 0, 1, rtoli, 0.0d+0) 1475 go to 700 1476 620 call xerrwv('lsode-- atol(i1) is r1 .lt. 0.0 ', 1477 1 40, 20, 1, 1, i, 0, 1, atoli, 0.0d+0) 1478 go to 700 1479 621 ewti = rwork(lewt+i-1) 1480 call xerrwv('lsode-- ewt(i1) (=r1) is .le. 0.0 ', 1481 1 40, 21, 1, 1, i, 0, 1, ewti, 0.0d+0) 1482 go to 700 1483 622 call xerrwv( 1484 1 'lsode-- tout (=r1) too close to t(=r2) ', 1485 1 60, 22, 1, 0, 0, 0, 2, tout, t) 1486 go to 700 1487 623 call xerrwv( 1488 1 'lsode-- itask (=i1) and tout (=r1) .gt. tcur - hu (= r2) ', 1489 1 60, 23, 1, 1, itask, 0, 2, tout, tp) 1490 go to 700 1491 624 call xerrwv( 1492 1 'lsode-- itask = 4 or 5 and tcrit (=r1) .gt. tcur (=r2) ', 1493 1 60, 24, 1, 0, 0, 0, 2, tcrit, tn) 1494 go to 700 1495 625 call xerrwv( 1496 1 'lsode-- itask = 4 or 5 and tcrit (=r1) .lt. tout (=r2)', 1497 1 60, 25, 1, 0, 0, 0, 2, tcrit, tout) 1498 go to 700 1499 626 call xerrwv('lsode-- initial precision required', 1500 1 50, 26, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1501 call xerrwv( 1502 1 'too high wrt machine precision tolsf (=r1)', 1503 1 60, 26, 1, 0, 0, 0, 1, tolsf, 0.0d+0) 1504 rwork(14) = tolsf 1505 go to 700 1506 627 call xerrwv('lsode-- problems in intdy. itask=i1,tout=r1', 1507 1 50, 27, 1, 1, itask, 0, 1, tout, 0.0d+0) 1508c 1509 700 if (illin .eq. 5) go to 710 1510 illin = illin + 1 1511 istate = -3 1512 return 1513 710 call xerrwv('lsode-- incorrect inputs', 1514 1 50, 302, 1, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1515c 1516 800 call xerrwv('lsode-- infinite loop ', 1517 1 50, 303, 2, 0, 0, 0, 0, 0.0d+0, 0.0d+0) 1518 return 1519c----------------------- end of subroutine lsode ----------------------- 1520 end 1521