1 Brief Description of 2 ODEPACK - A Systematized Collection of ODE Solvers 3 Double Precision Version 4 5 6 Alan C. Hindmarsh 7 Center for Applied Scientific Computing, L-561 8 Lawrence Livermore National Laboratory 9 Livermore, CA 94551, U.S.A. 10 11 20 June 2001 12 13 14 15Work performed under the auspices of the U.S. Department of Energy 16by the Lawrence Livermore National Laboratory under contract 17No. W-7405-Eng-48, and supported (formerly) by the DOE Office of 18Energy Research, Applied Mathematical Sciences Research Program. 19 20--------------------------------------------------------------------------- 21 22 ODEPACK is a collection of Fortran solvers for the initial value 23problem for ordinary differential equation systems. It consists of nine 24solvers, namely a basic solver called LSODE and eight variants of it -- 25LSODES, LSODA, LSODAR, LSODPK, LSODKR, LSODI, LSOIBT, and LSODIS. 26The collection is suitable for both stiff and nonstiff systems. It 27includes solvers for systems given in explicit form, dy/dt = f(t,y), 28and also solvers for systems given in linearly implicit form, 29A(t,y) dy/dt = g(t,y). Two of the solvers use general sparse matrix 30solvers for the linear systems that arise. Two others use iterative 31(preconditioned Krylov) methods instead of direct methods for these 32linear systems. The most recent addition is LSODIS, which solves 33implicit problems with general sparse treatment of all matrices involved. 34 35 The ODEPACK solvers are written in standard Fortran 77, with a few 36exceptions, and with minimal machine dependencies. There are separate 37double and single precision versions of ODEPACK. The actual solver 38names are those given above with a prefix of D- or S- for the double 39or single precision version, respectively, i.e. DLSODE/SLSODE, etc. 40Each solver consists of a main driver subroutine having the same name 41as the solver and some number of subordinate routines. For each 42solver, there is also a demonstration program, which solves one or two 43simple problems in a somewhat self-checking manner. 44 45 Recently, the ODEPACK solvers were upgraded to improve their 46portability in numerous ways. Among the improvements are (a) renaming 47of routines and Common blocks to distinguish double and single 48precision versions, (b) use of generic intrinsic function names, (c) 49elimination of the Block Data subprogram, (d) use of a portable 50routine to set the unit roundoff, and (e) passing of quoted strings to 51the error message handler. In addition, the prologue and internal 52comments were reformatted, and use mixed upper/lower case. Numerous 53minor corrections and improvements were also made. 54 55 The above upgrade operations were applied to LSODE earlier than they 56were to the rest of ODEPACK, and the two upgrades were done somewhat 57independently. As a result, some differences will be apparent in the 58source files of LSODE and the other solvers -- primarily in the 59formatting of the comment line prologue of the main driver routine. 60In Subroutines DLSODE/SLSODE and their subordinate routines, the 61prologue was written in "SLATEC format", while for the other solvers a 62more relaxed style was used. The differences are entirely cosmetic, 63however, and do not affect performance. 64 65 Documentation on the usage of each solver is provided in the 66initial block of comment lines in the source file, which (in most 67cases) includes a simple example. A demonstration program (in 68seperate double/single precision versions) is also available. 69 70 What follows is a summary of the capabilities of ODEPACK, comments 71about usage documentation, and notes about installing the collection. 72For additional documentation on ODEPACK, see also the papers [1], [2] 73(for LSODE), and [3] (for LSODPK and LSODKR), and in the references 74cited there. (However, the document [2] does not reflect the upgrade 75operations described above.) 76 77 78References: 79 80[1] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE Solvers," 81 in Scientific Computing, R. S. Stepleman et al. (eds.), North-Holland, 82 Amsterdam, 1983 (vol. 1 of IMACS Transactions on Scientific Computation), 83 pp. 55-64. 84 85[2] K. Radhakrishnan and A. C. Hindmarsh, "Description and Use of LSODE, 86 the Livermore Solver for Ordinary Differential Equations," LLNL 87 report UCRL-ID-113855, December 1993. 88 89[3] P. N. Brown and A. C. Hindmarsh, "Reduced Storage Matrix Methods 90 in Stiff ODE Systems," J. Appl. Math. & Comp., 31 (1989), pp.40-91. 91 92--------------------------------------------------------------------------- 93 94 I. Summary of the ODEPACK Solvers 95 96A. Solvers for explicitly given systems. 97 98For each of the following solvers, it is assumed that the ODEs are 99given explicitly, so that the system can be written in the form 100dy/dt = f(t,y), where y is the vector of dependent variables, and t is 101the independent variable. 102 1031. LSODE (Livermore Solver for Ordinary Differential Equations) is the 104 basic solver of the collection. It solves stiff and nonstiff systems 105 of the form dy/dt = f. In the stiff case, it treats the Jacobian 106 matrix df/dy as either a dense (full) or a banded matrix, and as either 107 user-supplied or internally approximated by difference quotients. 108 It uses Adams methods (predictor-corrector) in the nonstiff case, 109 and Backward Differentiation Formula (BDF) methods (the Gear methods) 110 in the stiff case. The linear systems that arise are solved by direct 111 methods (LU factor/solve). LSODE supersedes the older GEAR and GEARB 112 packages, and reflects a complete redesign of the user interface 113 and internal organization, with some algorithmic improvements. 114 1152. LSODES, written jointly with A. H. Sherman, solves systems dy/dt = f 116 and in the stiff case treats the Jacobian matrix in general sparse 117 form. It determines the sparsity structure on its own, or optionally 118 accepts this information from the user. It then uses parts of the 119 Yale Sparse Matrix Package (YSMP) to solve the linear systems that 120 arise, by a sparse (direct) LU factorization/backsolve method. 121 LSODES supersedes, and improves upon, the older GEARS package. 122 1233. LSODA, written jointly with L. R. Petzold, solves systems dy/dt = f 124 with a dense or banded Jacobian when the problem is stiff, but it 125 automatically selects between nonstiff (Adams) and stiff (BDF) 126 methods. It uses the nonstiff method initially, and dynamically 127 monitors data in order to decide which method to use. 128 1294. LSODAR, also written jointly with L. R. Petzold, is a variant of 130 LSODA with a rootfinding capability added. Thus it solves problems 131 dy/dt = f with dense or banded Jacobian and automatic method 132 selection, and at the same time, it finds the roots of any of a 133 set of given functions of the form g(t,y). This is often useful 134 for finding stop conditions, or for finding points at which a switch 135 is to be made in the function f. 136 1375. LSODPK, written jointly with Peter N. Brown, is a variant of LSODE 138 in which the direct solvers for the linear systems have been replaced 139 by a selection of four preconditioned Krylov (iterative) solvers. 140 The user must supply a pair of routine to evaluate, preprocess, and 141 solve the (left and/or right) preconditioner matrices. LSODPK also 142 includes an option for a user-supplied linear system solver to be used 143 without Krylov iteration. 144 1456. LSODKR is a variant of LSODPK with the addition of the same 146 rootfinding capability as in LSODAR, and also of automatic switching 147 between functional and Newton iteration. The nonlinear iteration 148 method-switching differs from the method-switching in LSODA and LSODAR, 149 but provides similar savings by using the cheaper method in the non-stiff 150 regions of the problem. LSODKR also improves on the Krylov methods in 151 LSODPK by offering the option to save and reuse the approximate Jacobian 152 data underlying the preconditioner. 153 154 155B. Solvers for linearly implicit systems. 156 157The following solvers treat systems in the linearly implicit form 158A(t,y) dy/dt = g(t,y), A = a square matrix, i.e. with the derivative 159dy/dt implicit, but linearly so. These solvers allow A to be 160singular, in which case the system is a differential-algebraic 161equation (DAE) system. In that case, the user must be very careful 162to supply a well-posed problem with consistent initial conditions. 163 1647. LSODI, written jointly with J. F. Painter, solves linearly implicit 165 systems in which the matrices involved (A, dg/dy, and d(A dy/dt)/dy) 166 are all assumed to be either dense or banded. LSODI supersedes the 167 older GEARIB solver and improves upon it in numerous ways. 168 1698. LSOIBT, written jointly with C. S. Kenney, solves linearly implicit 170 systems in which the matrices involved are all assumed to be 171 block-tridiagonal. Linear systems are solved by the LU method. 172 1739. LSODIS, written jointly with S. Balsdon, solves linearly implicit 174 systems in which the matrices involved are all assumed to be sparse. 175 Like LSODES, LSODIS either determines the sparsity structure or 176 accepts it from the user, and uses parts of the Yale Sparse Matrix 177 Package to solve the linear systems that arise, by a direct method. 178 179--------------------------------------------------------------------------- 180 181 II. Usage Documentation 182 183 Each of the solvers in the ODEPACK collection is headed by a 184user-callable driver subroutine, with the same name as the solver 185(DLSODE, etc.). The call sequence of the driver routine includes the 186names of one or more user-supplied subroutines that define the ODE 187system, and various other problem and solution parameters. Complete 188user documentation is given in the initial block of comment lines 189(the prologue) in the driver routine. In each case, this prologue is 190organized as follows: 191 192 * Summary of Usage (short, for standard modes of use) 193 * Example Problem, with code and output (except for LSODPK and LSODKR) 194 * Full Description of User Interface, further divided as follows: 195 1. Call sequence description (including optional inputs/outputs) 196 2. Optionally callable routines 197 3. Descriptions of internal Common blocks 198 4. Optionally user-replaceable routines 199 * Revision History, showing date written and dates of revisions 200 * Other Routines, a list of all subordinate routines for the solver 201 202 First-time users should read only the Summary of Usage and look at the 203the Example Problem (or demonstration program), then later refer to the 204Full Description if and when more details or nonstandard options are needed. 205 206--------------------------------------------------------------------------- 207 208 III. Installation Notes 209 2101. In addition to this document, the double precision version of ODEPACK 211consists of three source files, plus a demonstration program file. 212The solver source files are organized as follows: 213 214 opkdmain.f = Main Source File, consisting of the driver subroutine 215 for each of the nine solvers 216 217 opkda1.f = First Auxiliary Source File, consisting of subordinate 218 routines for the solvers 219 220 opkda2.f = Second Auxiliary Source File, consisting of subordinate 221 routines which may already reside on the user's system 222 (See Notes 2 and 3 below.) 223 224The demonstration program file is: 225 226 opkddemos = a merge of the nine demonstration programs, with the 227 source for each followed by a sample output 228 2292. The Second Auxiliary Source File includes the routines from the 230LINPACK and BLAS collections that are needed by the solvers (and by two 231of the demonstration programs), for the solution of dense and banded 232linear systems and associated basic linear algebra operations. 233These routine are: 234 From LINPACK: DGEFA, DGESL, DGBFA, DGBSL 235 From the BLAS: DAXPY, DCOPY, DDOT, DSCAL, DNRM2, IDAMAX 236If your computer system already has these routines, and especially if it 237has machine-optimized versions, the copies provided here can be discarded. 238 2393. The Second Auxiliary Source File includes a set of five routines -- 240XERRWD, XSETUN, XSETF, IXSAV, IUMACH -- which handle error messages 241from the solvers. This set is in fact a reduced version (sufficient 242for the needs of ODEPACK) of a much larger error handling package from 243the SLATEC Library. If your computer system already has the full 244SLATEC error handler, the version provided here can be discarded. If 245the reduced version is used, its machine-dependent features should be 246checked first; see comments in Subroutine XERRWD. 247 2484. ODEPACK contains a few instances where ANSI Fortran 77 is violated: 249 (a) In various places in the LSODES and LSODIS solvers, a call to a 250 subroutine has a subscripted real array as an argument where the 251 subroutine called expects an integer array. Calls of this form 252 occur in Subroutine DLSODES (to DSTODE), in DIPREP (to DPREP), 253 in Subroutine DLSODIS (to DSTODI), and in DIPREPI (to DPREPI). 254 Another such call occurs in the DLSODES demonstration program, 255 from the main program to Subroutine SSOUT. This is done in order 256 to use work space in an efficient manner, as the same space is 257 sometimes used for real work space and sometimes for integer work 258 space. If your compiler does not accept this feature, one possible 259 way to get the desired result is to compile the called routines 260 and calling routines in separate jobs, and then combine the binary 261 modules in an appropriate manner. If this procedure is still not 262 acceptable under your system, it will be necessary to radically 263 alter the structure of the array RWORK within the LSODES or LSODIS 264 solver package. (See also Note 5 below.) 265 (b) Each ODEPACK solver treats the arguments NEQ, Y, RTOL, and ATOL 266 as arrays, even though the length may be only 1. Moreover, 267 except for Y, the usage instructions say that these arguments 268 may be either arrays or scalars. If your system does not allow 269 such a mismatch, then the documentation of these arguments 270 should be changed accordingly. 271 2725. For maximum storage economy, the LSODES and LSODIS solvers make use 273of the real to integer wordlength ratio. This is assumed to be an 274integer L such that if a real array R and an integer array M occupy 275the same space in memory, R(1) having the same bit address as M(1), 276then R(I) has the same address as M((I-1)*L+1). This ratio L is 277usually 2 for double precision, and this is the value used in the 278double precision version supplied. If this value is incorrect, it 279needs to be changed in two places: 280 (a) The integer LENRAT is DATA-loaded in Subroutines DLSODES and 281 DLSODIS to this ratio, shortly below the prologue. 282 (b) The integer LRATIO is DATA-loaded in Subroutine CDRV to this 283 ratio, shortly below the prologue of that routine. 284(See comments in both places.) If the ratio is not an integer, use 285the greatest integer not exceeding the ratio. 286 2876. For installation of ODEPACK on a Cray computer, the source files 288supplied include compiler directives for the CFT compiler. These have 289the form CDIR$ IVDEP and occur prior to certain loops that involve 290subscript shifts (and would otherwise not be vectorized). These 291directives are (or should be) treated as comments by any other compiler. 292 2937. On first obtaining ODEPACK, the demonstration programs should be 294compiled and executed prior to any other use of the solvers. In most 295cases, these excercise all of the major method options in each solver, 296and are self-checking. (In the case of LSODPK and LSODKR, the 297demonstration programs are not self-checking, and for LSODKR only one 298major method option is used.) In any case, the output can be compared 299with the sample output supplied, which was generated from the double 300precision version of ODEPACK on a 32-bit computer. When comparing your 301output with that supplied, differences of 10-20% in the final values of 302the various statistical counters can be attributed to differences in 303the roundoff properties of different computer systems. 304 3058. If some subset of the whole ODEPACK collection is desired, without 306unneeded routines, the appropriate routines must be extracted 307accordingly. The following lists give the routines needed for the 308double precision version of each solver. 309 310 The DLSODE solver consists of the routines 311 DLSODE, DINTDY, DSTODE, DCFODE, DPREPJ, DSOLSY, DEWSET, DVNORM, DSRCOM, 312 DGEFA, DGESL, DGBFA, DGBSL, DAXPY, DSCAL, DDOT, IDAMAX, 313 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 314 315 The DLSODES solver consists of the routines 316 DLSODES, DIPREP, DPREP, JGROUP, ADJLR, CNTNZU, DINTDY, DSTODE, DCFODE, 317 DPRJS, DSOLSS, DEWSET, DVNORM, DSRCMS, 318 ODRV, MD, MDI, MDM, MDP, MDU, SRO, 319 CDRV, NROC, NSFC, NNFC, NNSC, NNTC, 320 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 321 322 The DLSODA solver consists of the routines 323 DLSODA, DINTDY, DSTODA, DCFODE, DPRJA, DSOLSY, DEWSET, 324 DMNORM, DFNORM, DBNORM, DSRCMA, 325 DGEFA, DGESL, DGBFA, DGBSL, DAXPY, DSCAL, DDOT, IDAMAX, 326 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 327 328 The DLSODAR solver consists of the routines 329 DLSODAR, DRCHEK, DROOTS, DINTDY, DSTODA, DCFODE, DPRJA, DSOLSY, DEWSET, 330 DMNORM, DFNORM, DBNORM, DSRCAR, 331 DGEFA, DGESL, DGBFA, DGBSL, DAXPY, DSCAL, DDOT, DCOPY, IDAMAX, 332 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 333 334 The DLSODPK solver consists of the routines 335 DLSODPK, DINTDY, DEWSET, DVNORM, DSTODPK, DCFODE, DPKSET, DSOLPK, 336 DSPIOM, DATV, DORTHOG, DHEFA, DHESL, DSPIGMR, DHEQR, DHELS, 337 DPCG, DPCGS, DATP, DUSOL, DSRCPK, 338 DAXPY, DSCAL, DCOPY, DDOT, DNRM2, IDAMAX, 339 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 340 341 The DLSODKR solver consists of the routines 342 DLSODKR, DRCHEK, DROOTS, DLHIN, DINTDY, DEWSET, DVNORM, DSTOKA, 343 DCFODE, DSETPK, DSOLPK, DSPIOM, DATV, DORTHOG, DHEFA, DHESL, DSPIGMR, 344 DHEQR, DHELS, DPCG, DPCGS, DATP, DUSOL, DSRCKR, 345 DAXPY, DSCAL, DCOPY, DDOT, DNRM2, IDAMAX, 346 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 347 348 The DLSODI solver consists of the routines 349 DLSODI, DAINVG, DINTDY, DSTODI, DCFODE, DPREPJI, DSOLSY, DEWSET, 350 DVNORM, DSRCOM, 351 DGEFA, DGESL, DGBFA, DGBSL, DAXPY, DSCAL, DDOT, IDAMAX, 352 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 353 354 The DLSOIBT solver consists of the routines 355 DLSOIBT, DAIGBT, DINTDY, DSTODI, DCFODE, DPJIBT, DSLSBT, DEWSET, 356 DVNORM, DSRCOM, DDECBT, DSOLBT, 357 DGEFA, DGESL, DAXPY, DSCAL, DDOT, IDAMAX, 358 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 359 360 The DLSODIS solver consists of the routines 361 DLSODIS, DAINVGS, DIPREPI, DPREPI, JGROUP, ADJLR, CNTNZU, DINTDY, 362 DSTODI, DCFODE, DPRJIS, DSOLSS, DEWSET, DVNORM, DSRCMS, 363 ODRV, MD, MDI, MDM, MDP, MDU, SRO, 364 CDRV, NROC, NSFC, NNFC, NNSC, NNTC, 365 DUMACH, XERRWD, XSETUN, XSETF, IXSAV, IUMACH 366