Brief Description of
ODEPACK - A Systematized Collection of ODE Solvers
Alan C. Hindmarsh
Computing & Mathematics Research Division, L-316
Lawrence Livermore National Laboratory
Livermore, CA 94550, U.S.A.
19 November 1985
Work performed under the auspices of the U.S. Department of Energy
by the Lawrence Livermore National Laboratory under contract
No. W-7405-Eng-48, and supported by the DOE Office of Energy Research,
Applied Mathematical Sciences Research Program.
--------------------------------------------------------------------------------
ODEPACK is a collection of Fortran solvers for the initial value
problem for ordinary differential equation (ODE) systems. It currently
includes six solvers, suitable for both stiff and nonstiff systems,
and includes solvers for systems given in linearly implicit form as well
as solvers for systems given in explicit form.
The solvers are written in Fortran IV (1966 ANSI Fortran) with a few
exceptions, and with minimal machine dependencies. Each solver consists of
a driver having the same name as the solver (i.e. LSODE etc.), and some
number of subordinate routines.
What follows is a summary of the capabilities of ODEPACK, comments
about usage documentation, and notes about installing the collection.
Further details are available in the reference given at the end, and in the
references cited there.
--------------------------------------------------------------------------------
I. Summary of the ODEPACK Solvers
A. Solvers for explicitly given systems.
In the solvers below, it is assumed that the ODE's are given
explicitly, so that the system can be written in the form
dy/dt = f(t,y) ,
where y is the vector of dependent variables, and t is the independent
variable.
1. LSODE (Livermore Solver for Ordinary Differential Equations) is the
basic solver of the collection. It solves stiff and nonstiff systems
of the form dy/dt = f. In the stiff case, it treats the Jacobian matrix
df/dy as either a full or a banded matrix, and as either user-supplied
or internally approximated by difference quotients. It uses Adams methods
(predictor-corrector) in the nonstiff case, and Backward Differentiation
Formula (BDF) methods in the stiff case. The linear systems that arise
are solved by direct methods (LU factor/solve). LSODE supersedes the older
GEAR and GEARB packages, and reflects a complete redesign of the user
interface and internal organization, with some algorithmic improvements.
2. LSODES, written jointly with A. H. Sherman, solves systems dy/dt = f
and in the stiff case treats the Jacobian matrix in general sparse
form. It determines the sparsity structure on its own (or optionally
accepts this information from the user), and uses parts of the Yale Sparse
Matrix Package (YSMP) to solve the linear systems that arise.
LSODES supersedes, and improves upon, the older GEARS package.
3. LSODA, written jointly with L. R. Petzold, solves systems dy/dt = f
with a full or banded Jacobian when the problem is stiff, but it
automatically selects between nonstiff (Adams) and stiff (BDF) methods.
It uses the nonstiff method initially, and dynamically monitors data
in order to decide which method to use.
4. LSODAR, also written jointly with L. R. Petzold, is a variant of LSODA
with a rootfinding capability added. Thus it solves problems dy/dt = f
with full or banded Jacobian and automatic method selection, and at
the same time, it finds the roots of any of a set of given functions
of the form g(t,y). This is often useful for finding stop conditions
or points at which switches are to be made in the function f.
B. Solvers for linearly implicit systems.
In the solvers below, it is assumed that the derivative dy/dt is
implicit, but occurs linearly, so that the system can be written
A(t,y) dy/dt = g(t,y) ,
where A is a square matrix. These solvers allow A to be singular,
in which case the system is a differential-algebraic system, but in that
case users must be very careful to supply a well-posed problem with
consistent initial conditions.
5. LSODI, written jointly with J. F. Painter, solves linearly implicit
systems in which the matrices involved (A, dg/dy, and d(A dy/dt)/dy) are
all assumed to be either full or banded. LSODI supersedes the older
GEARIB solver and improves upon it in numerous ways.
6. LSOIBT, written jointly with C. S. Kenney, solves linearly implicit
systems in which the matrices involved are all assumed to be
block-tridiagonal. Linear systems are solved by the LU method.
--------------------------------------------------------------------------------
II. Usage Documentation
Each of the solvers in the ODEPACK collection is headed by a
user-callable driver subroutine, whose call sequence includes the names
of one or more user-supplied subroutines that define the ODE system.
Complete user documentation is given in the initial block of comment
cards (the prologue) in the driver routine. In each case, this prologue
is organized as follows:
Summary of Usage (short, for standard modes of use)
Example Problem (with code and output)
Full Description of User Interface, further divided as follows:
I. Call sequence description (including optional inputs/outputs)
II. Optionally callable routines
III. Descriptions of internal Common blocks
IV. Optionally user-replaceable routines
First-time users should read only the Summary of Usage and look
at the Example Problem, then later refer to the Full Description if and
when more details or nonstandard options are needed.
Depending on the user environment, it may be desirable to create
separate copies of these prologues as user manuals. (The installation notes
that follow should be checked first, however, to see if usage instructions
are affected by changes made upon istallation.)
--------------------------------------------------------------------------------
III. Installation Notes
1. ODEPACK is being made available in separate single and double
precision versions. For each precision, a file is provided containing
appropriate subroutines and function routines in source form, along
with a file containing demonstration (quick check) programs.
It is expected that any one user will be interested in only one of
the two precisions. Therefore, not all of the routine names in the two
precisions were made distinct. If instead both precisions of one or more
of the solvers are to be combined as a library, then renaming must be done
for all routines which involve real numbers and which have non-unique
names initially (e.g. SLSODE and DLSODE instead of LSODE).
2. These source files are complete except for needed routines from the
LINPACK and BLAS collections and machine constant routines. These are:
From LINPACK: SGEFA, SGESL, SGBFA, SGBSL (in single precision solvers)
DGEFA, DGESL, DGBFA, DGBSL (in double precision solvers)
From the BLAS: SAXPY, SCOPY, SDOT, SSCAL, ISAMAX (single precision)
DAXPY, DCOPY, DDOT, DSCAL, IDAMAX (double precision)
Machine constant routines: R1MACH (in S.P.), D1MACH (in D.P.)
(used only to provide the unit roundoff)
3. The source files (for each precision) include a small error handling
package XERRWV/XSETUN/XSETF. This is a reduced version of the much larger
SLATEC Error Handling Package, and is sufficient for the needs of ODEPACK.
If the SLATEC version is available, the reduced version can be discarded.
However, in that case it may be necessary to first execute a call
CALL XSETF(1)
before calling any ODEPACK solver, in order to insure that recoverable
errors are not treated as fatal. Alternatively, set the default value of
IFLAG to 1 in the SLATEC error package. If the reduced version is used,
its machine-dependent features should be checked first (see comments in
Subroutine XERRWV).
4. In each solver, there are four integer variables, in two internal labeled
Common blocks, which need to be loaded with DATA statements. (They
can vary during execution, and are in Common to assure their
retention between calls.) This is legal in ANSI Fortran only if
done in a Block Data subprogram, and the double precision source
file has a Block Data subprogram for this purpose. However, because
Block Data subprograms can be difficult to install in libraries,
and because many compilers allow such DATA statements in subroutines,
the single precision source file is being supplied with
two such DATA statements. If, for your system, the way in which
this is handled is incorrect in the version you are installing,
it is easy to change to the other way. The locations of these
DATA statements in a version without Block Data are just after the
initial type and Common declarations in each driver subroutine and
in XERRWV. In each driver, ILLIN and NTREP are DATA-loaded as 0.
In XERRWV, MESFLG is loaded as 1 and LUNIT is loaded as the
appropriate default logical unit number. The Block Data subprogram
(if used) must of course have copies of the Common declarations,
and also a type declaration in the double precision case. If
necessary for uniqueness, name the subprogram.
5. ODEPACK contains a few instances where ANSI Fortran (1966) is violated:
(a) Subscripts of the form I + J, I - J, I + J + const,
and I + J - const occur in various routines.
In the sparse matrix routines, there are subscripted subscripts,
i.e. subsripts of the form I(J), I(J)+1, -I(J), I(J+K), and
(in NNFC only) I(J(K)).
If any of these forms is unacceptable to your compiler, make the
obvious changes to the source code accordingly. (But this should
be done only where necessary, to avoid loss of efficiency.)
(b) The intrinsic function DFLOAT (conversion from integer to double
precision) is used by ODEPACK, but is not required in ANSI Fortran.
Thus it may have to be supplied separately for the double precision
versions.
(c) In various places in the LSODES solver, a call to a
subroutine has a subscripted real array as an argument where
the subroutine called has an integer array. Calls of this form
occur in Subroutine LSODES (to STODE) and in IPREP (to PREP).
Another such call occurs in the LSODES demonstration program, from
the main program to Subroutine SSOUT. This is done in order
to use work space in an efficient manner, as the same space
is sometimes used for real work space and sometimes for integer
work space. If your compiler does not accept this feature,
one possible way to get the desired result is to compile the
called routines and calling routines in separate jobs, and then
combine the binary modules in an appropriate manner. If this
procedure is still not acceptable under your system, it will
be necessary to alter radically the structure of the array RWORK
within the LSODES package. (See also the note below about LSODES.)
(d) Each ODEPACK solver treats the arguments NEQ, Y, RTOL, and ATOL as
arrays, even though the length may be only 1. Moreover, except for Y,
the usage instructions say that these arguments may be either arrays
or scalars. If your system does not allow such a mismatch, then the
documentation of these arguments should be changed accordingly.
6. For maximum storage economy, the LSODES solver makes
use of the real to integer wordlength ratio. This is assumed to be
an integer L such that if a real array R and an integer array M
occupy the same space in memory, R(1) having the same bit address
as M(1), then R(I) has the same address as M((I-1)*L+1).
This ratio L is usually 1 for single precision and 2 for double
precision, and these are the values used in the single and double
precision versions supplied, respectively. If the value supplied is
incorrect, it needs to be changed in two places:
(a) The integer LENRAT is DATA-loaded in Subroutine LSODES
to this ratio, shortly below the prologue.
(b) The integer LRATIO is DATA-loaded in Subroutine CDRV to this
ratio, shortly below the prologue of that routine.
(See comments in both places.) If the ratio is not an integer, use the
greatest integer not exceeding the ratio.
7. For installation of ODEPACK on a Cray computer, the source files
supplied include compiler directives for the CFT compiler.
These have the form CDIR$ IVDEP and occur prior to certain loops
that involve subscript shifts (and would otherwise not be vectorized).
These directives are (or should be) treated as comments by any
other compiler.
8. On first obtaining ODEPACK, the demonstration programs should be
compiled and executed prior to any other use of the solvers. These
excercise all of the major method options in each solver, and are
self-checking.
9. If some subset of the whole ODEPACK collection is desired, without
unneeded routines, the appropriate routines must be extracted accordingly.
The following lists give the routines needed for each solver.
The LSODE solver consists of the routines
LSODE, INTDY, STODE, CFODE, PREPJ, SOLSY, EWSET, VNORM, SVCOM, RSCOM,
SGEFA or DGEFA, SGESL or DGESL, SGBFA or DGBFA, SGBSL or DGBSL,
SAXPY or DAXPY, SSCAL or DSCAL, ISAMAX or IDAMAX, SDOT or DDOT,
R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only)
The LSODES solver consists of the routines
LSODES, IPREP, PREP, JGROUP, ADJLR, CNTNZU, INTDY, STODE, CFODE,
PRJS, SLSS, EWSET, VNORM, SVCMS, RSCMS,
ODRV, MD, MDI, MDM, MDP, MDU, SRO,
CDRV, NROC, NSFC, NNFC, NNSC, NNTC,
R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only)
The LSODA solver consists of the routines
LSODA, INTDY, STODA, CFODE, PRJA, SOLSY, EWSET,
VMNORM, FNORM, BNORM, SVCMA, RSCMA,
SGEFA or DGEFA, SGESL or DGESL, SGBFA or DGBFA, SGBSL or DGBSL,
SAXPY or DAXPY, SSCAL or DSCAL, ISAMAX or IDAMAX, SDOT or DDOT,
R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only)
The LSODAR solver consists of the routines
LSODAR, RCHEK, ROOTS, INTDY, STODA, CFODE, PRJA, SOLSY, EWSET,
VMNORM, FNORM, BNORM, SVCAR, RSCAR,
SGEFA or DGEFA, SGESL or DGESL, SGBFA or DGBFA, SGBSL or DGBSL,
SAXPY or DAXPY, SSCAL or DSCAL, ISAMAX or IDAMAX, SDOT or DDOT,
SCOPY or DCOPY, R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only)
The LSODI solver consists of the routines
LSODI, AINVG, INTDY, STODI, CFODE, PREPJI, SOLSY, EWSET, VNORM,
SVCOM, RSCOM,
SGEFA or DGEFA, SGESL or DGESL, SGBFA or DGBFA, SGBSL or DGBSL,
SAXPY or DAXPY, SSCAL or DSCAL, ISAMAX or IDAMAX, SDOT or DDOT,
R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only).
The LSOIBT solver consists of the routines
LSOIBT, AIGBT, INTDY, STODI, CFODE, PJIBT, SLSBT, EWSET, VNORM,
SVCOM, RSCOM,
DECBT, SOLBT, SGEFA or DGEFA, SGESL or DGESL,
SAXPY or DAXPY, SSCAL or DSCAL, ISAMAX or IDAMAX, SDOT or DDOT,
R1MACH or D1MACH, XERRWV, XSETUN, XSETF,
and a Block Data subprogram (double precision version only).
--------------------------------------------------------------------------------
Reference
A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE Solvers,"
in Scientific Computing, R. S. Stepleman et al. (eds.), North-Holland,
Amsterdam, 1983 (vol. 1 of IMACS Transactions on Scientific Computation),
pp. 55-64. Also available as LLNL Report UCRL-88007, August 1982.
--------------------------------------------------------------------------------