dvode: variable-coefficient ordinary differential equation solver, with fixed-leading-coefficient implementation.
dvode solves the initial value problem for stiff or nonstiff systems of first order odes:
or, in component form:
dvode is a package based on the episode and episodeb packages, and
on the odepack user interface standard, with minor modifications.
if the integration of a given problem by dvode is to be interrrupted and then later continued, such as when restarting an interrupted run or alternating between two or more ode problems, the user should save, following the return from the last dvode call prior to the interruption, the contents of the call sequence variables and internal variables, and later restore these values before the next dvode call for that problem. to save and restore the variables, use subroutine dvsrco.
in addition, if non-default values for either lun or mflag are
desired, an extra call to xsetun and/or xsetf should be made just
before continuing the integration.
Note
the legality of input parameters will be thoroughly checked on the
initial call for the problem, but not checked thereafter unless a
change in input parameters is flagged by istate = 3 in the input.
Note
the work arrays must not be altered between calls to dvode
for the same problem, except possibly for the conditional and
optional input, and except for the last 3*neq words of rwork.
the latter space is used for internal scratch space, and so is
available for use by the user outside dvode between calls, if
desired (but not for use by f or jac).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dvode_t), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | neq |
the size of the ode system (number of first order
ordinary differential equations).
|
||
| real(kind=wp), | intent(inout) | :: | y(*) |
a real array for the vector of dependent variables, of
length this array is passed as the |
||
| real(kind=wp), | intent(inout) | :: | t |
the independent variable. in the input, |
||
| real(kind=wp), | intent(in) | :: | tout |
the next value of when starting the problem ( if if |
||
| integer, | intent(in) | :: | itol |
an indicator for the type of error control.
1 or 2 according as |
||
| real(kind=wp), | intent(in) | :: | rtol(*) |
a relative error tolerance parameter, either a scalar or
an array of length |
||
| real(kind=wp), | intent(in) | :: | atol(*) |
an absolute error tolerance parameter, either a scalar or
an array of length the input parameters and the rms-norm (root-mean-square norm) here is
when either of these parameters is a scalar, it need not be dimensioned in the user's calling program. if none of the above choices (with if global errors are to be estimated by making a repeated
run on the same problem with smaller tolerances, then all
components of use |
||
| integer, | intent(in) | :: | itask |
an index specifying the task to be performed. input only. itask has the following values and meanings:
note: if itask = 4 or 5 and the solver reaches tcrit (within roundoff), it will return t = tcrit (exactly) to indicate this (unless itask = 4 and tout comes before tcrit, in which case answers at t = tout are returned first). |
||
| integer, | intent(inout) | :: | istate |
an index used for input and output to specify the the state of the calculation. in the input, the values of istate are as follows:
note: a preliminary call with tout = t is not counted as a first call here, as no initialization or checking of input is done. (such a call is sometimes useful to include the initial conditions in the output.) thus the first call for which tout /= t requires istate = 1 in the input. in the output, istate has the following values and meanings:
note: since the normal output value of istate is 2, it does not need to be reset for normal continuation. also, since a negative input value of istate will be regarded as illegal, a negative output value requires the user to change it, and possibly other input, before calling the solver again. |
||
| integer, | intent(in) | :: | iopt |
an integer flag to specify whether or not any optional input is being used on this call. the optional input is listed separately below.
|
||
| real(kind=wp) | :: | rwork(lrw) |
a real working array.
the length of
(see the
the first 20 words of
the following word in
The following optional input requires
the following optional outputs
below are quantities related to the performance of dvode
which are available to the user.
except where stated otherwise, all of this output is defined
on any successful return from dvode, and on any return with
istate = -1, -2, -4, -5, or -6. on an illegal input return
(istate = -3), they will be unchanged from their existing values
(if any), except possibly for
the following two arrays are segments of the
|
|||
| integer, | intent(in) | :: | lrw |
the length of the array rwork, as declared by the user. (this will be checked by the solver.) |
||
| integer | :: | iwork(liw) |
an integer work array. the length of iwork must be at least:
the first 30 words of the following 2 words in
these are the lower and upper
half-bandwidths, respectively, of the
banded jacobian, excluding the main diagonal.
the band is defined by the matrix locations
The following optional input requires
as optional additional output from dvode, the variables listed
below are quantities related to the performance of dvode
which are available to the user.
except where stated otherwise, all of this output is defined
on any successful return from dvode, and on any return with
istate = -1, -2, -4, -5, or -6. on an illegal input return
(istate = -3), they will be unchanged from their existing values
(if any), except possibly for
|
|||
| integer, | intent(in) | :: | liw |
the length of the array iwork, as declared by the user. (this will be checked by the solver.) |
||
| integer, | intent(in) | :: | mf |
method flag. standard values are:
the complete set of legal values of mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, -11, -12, -14, -15, -21, -22, -24, -25.
if miter = 1 or 4, the user must supply a subroutine |
subroutine dvode(me,neq,y,t,tout,itol,rtol,atol,itask,istate,iopt, & rwork,lrw,iwork,liw,mf) class(dvode_t),intent(inout) :: me real(wp),intent(inout) :: y(*) !! a real array for the vector of dependent variables, of !! length `neq` or more. used for both input and output on the !! first call (`istate = 1`), and only for output on other calls. !! on the first call, `y` must contain the vector of initial !! values. in the output, `y` contains the computed solution !! evaluated at `t`. if desired, the `y` array may be used !! for other purposes between calls to the solver. !! !! this array is passed as the `y` argument in all calls to !! `f` and `jac`. real(wp),intent(inout) :: t !! the independent variable. in the input, `t` is used only on !! the first call, as the initial point of the integration. !! in the output, after each call, `t` is the value at which a !! computed solution `y` is evaluated (usually the same as `tout`). !! on an error return, `t` is the farthest point reached. real(wp),intent(in) :: tout !! the next value of `t` at which a computed solution is desired. !! !! when starting the problem (`istate = 1`), `tout` may be equal !! to `t` for one call, then should `/= t` for the next call. !! for the initial `t`, an input value of `tout /= t` is used !! in order to determine the direction of the integration !! (i.e. the algebraic sign of the step sizes) and the rough !! scale of the problem. integration in either direction !! (forward or backward in `t`) is permitted. !! !! if `itask = 2 or 5` (one-step modes), `tout` is ignored after !! the first call (i.e. the first call with `tout /= t`). !! otherwise, `tout` is required on every call. !! !! if `itask = 1, 3, or 4`, the values of `tout` need not be !! monotone, but a value of `tout` which backs up is limited !! to the current internal `t` interval, whose endpoints are !! `tcur - hu` and `tcur`. (see optional output, below, for !! `tcur` and `hu`.) real(wp),intent(in) :: rtol(*) !! a relative error tolerance parameter, either a scalar or !! an array of length `neq`. see description under `atol`. real(wp),intent(in) :: atol(*) !! an absolute error tolerance parameter, either a scalar or !! an array of length `neq`. !! !! the input parameters `itol`, `rtol`, and `atol` determine !! the error control performed by the solver. the solver will !! control the vector `e = (e(i))` of estimated local errors !! in `y`, according to an inequality of the form !!``` !! rms-norm of ( e(i)/ewt(i) ) <= 1, !! where ewt(i) = rtol(i)*abs(y(i)) + atol(i), !!``` !! and the rms-norm (root-mean-square norm) here is !! `rms-norm(v) = sqrt(sum v(i)**2 / neq)`. here `ewt = (ewt(i))` !! is a vector of weights which must always be positive, and !! the values of `rtol` and `atol` should all be non-negative. !! the following table gives the types (scalar/array) of !! `rtol` and `atol`, and the corresponding form of `ewt(i)`. !! !!``` !! itol rtol atol ewt(i) !! 1 scalar scalar rtol*abs(y(i)) + atol !! 2 scalar array rtol*abs(y(i)) + atol(i) !! 3 array scalar rtol(i)*abs(y(i)) + atol !! 4 array array rtol(i)*abs(y(i)) + atol(i) !!``` !! !! when either of these parameters is a scalar, it need not !! be dimensioned in the user's calling program. !! !! if none of the above choices (with `itol`, `rtol`, and `atol` !! fixed throughout the problem) is suitable, more general !! error controls can be obtained by substituting !! user-supplied routines for the setting of `ewt` and/or for !! the norm calculation. !! !! if global errors are to be estimated by making a repeated !! run on the same problem with smaller tolerances, then all !! components of `rtol` and `atol` (i.e. of `ewt`) should be scaled !! down uniformly. !! !! use `rtol = 0.0` for pure absolute error control, and !! use `atol = 0.0` (or `atol(i) = 0.0`) for pure relative error !! control. caution: actual (global) errors may exceed these !! local tolerances, so choose them conservatively. integer,intent(in) :: lrw !! the length of the array rwork, as declared by the user. !! (this will be checked by the solver.) real(wp) :: rwork(lrw) !! a real working array. !! the length of `rwork` must be at least !! `20 + nyh*(maxord + 1) + 3*neq + lwm` where: !! !! * `nyh` = the initial value of neq, !! * `maxord` = 12 (if `meth` = 1) or 5 (if `meth` = 2) (unless a !! smaller value is given as an optional input), !! * `lwm` = length of work space for matrix-related data: !! !! * `lwm = 0 ` if `miter` = 0, !! * `lwm = 2*neq**2 + 2 ` if `miter` = 1 or 2, and `mf`>0, !! * `lwm = neq**2 + 2 ` if `miter` = 1 or 2, and `mf`<0, !! * `lwm = neq + 2 ` if `miter` = 3, !! * `lwm = (3*ml+2*mu+2)*neq + 2` if `miter` = 4 or 5, and `mf`>0, !! * `lwm = (2*ml+mu+1)*neq + 2 ` if `miter` = 4 or 5, and `mf`<0. !! !! (see the `mf` description for `meth` and `miter`.) !! thus if maxord has its default value and `neq` is constant, !! this length is: !! !! * `20 + 16*neq ` for `mf` = 10, !! * `22 + 16*neq + 2*neq**2 ` for `mf` = 11 or 12, !! * `22 + 16*neq + neq**2 ` for `mf` = -11 or -12, !! * `22 + 17*neq ` for `mf` = 13, !! * `22 + 18*neq + (3*ml+2*mu)*neq` for `mf` = 14 or 15, !! * `22 + 17*neq + (2*ml+mu)*neq ` for `mf` = -14 or -15, !! * `20 + 9*neq ` for `mf` = 20, !! * `22 + 9*neq + 2*neq**2 ` for `mf` = 21 or 22, !! * `22 + 9*neq + neq**2 ` for `mf` = -21 or -22, !! * `22 + 10*neq ` for `mf` = 23, !! * `22 + 11*neq + (3*ml+2*mu)*neq` for `mf` = 24 or 25. !! * `22 + 10*neq + (2*ml+mu)*neq ` for `mf` = -24 or -25. !! !! the first 20 words of `rwork` are reserved for conditional !! and optional input and optional output. !! !! `rwork` can also used for conditional and !! optional input and optional output. !! !! the following word in `rwork` is a conditional input: !! !! * `rwork(1) = tcrit` = critical value of `t` which the solver !! is not to overshoot. required if `itask` is !! 4 or 5, and ignored otherwise. (see `itask`.) !! !! The following optional input requires `iopt = 1`, and in that !! case all of this input is examined. a value of zero for any !! of these optional input variables will cause the default value to be !! used. thus to use a subset of the optional input, simply preload !! locations 5 to 10 in `rwork` to 0.0, and !! then set those of interest to nonzero values: !! !! * `rwork(5) = h0` = the step size to be attempted on the first step. !! the default value is determined by the solver. !! !! * `rwork(6) = hmax` = the maximum absolute step size allowed. !! the default value is infinite. !! !! * `rwork(7) = hmin` = the minimum absolute step size allowed. !! the default value is 0. (this lower bound is not !! enforced on the final step before reaching `tcrit` !! when `itask = 4 or 5`.) !! !! the following optional outputs !! below are quantities related to the performance of [[dvode]] !! which are available to the user. !! except where stated otherwise, all of this output is defined !! on any successful return from [[dvode]], and on any return with !! istate = -1, -2, -4, -5, or -6. on an illegal input return !! (istate = -3), they will be unchanged from their existing values !! (if any), except possibly for `tolsf`. !! on any error return, output relevant to the error will be defined, !! as noted below: !! !! * `rwork(11) = hu` = the step size in t last used (successfully). !! !! * `rwork(12) = hcur` = the step size to be attempted on the next step. !! !! * `rwork(13) = tcur` = the current value of the independent variable !! which the solver has actually reached, i.e. the !! current internal mesh point in t. in the output, !! tcur will always be at least as far from the !! initial value of t as the current argument t, !! but may be farther (if interpolation was done). !! !! * `rwork(14) = tolsf` = a tolerance scale factor, greater than 1.0, !! computed when a request for too much accuracy was !! detected (istate = -3 if detected at the start of !! the problem, istate = -2 otherwise). if itol is !! left unaltered but rtol and atol are uniformly !! scaled up by a factor of tolsf for the next call, !! then the solver is deemed likely to succeed. !! (the user may also ignore tolsf and alter the !! tolerance parameters in any other way appropriate.) !! !! the following two arrays are segments of the `rwork` array which !! may also be of interest to the user as optional output: !! !! * `rwork(21:) = yh` = the nordsieck history array, of size `nyh` by !! `(nqcur + 1)`, where `nyh` is the initial value !! of `neq`. for `j = 0,1,...,nqcur`, column `j+1` !! of `yh` contains `hcur**j/factorial(j)` times !! the `j-th` derivative of the interpolating !! polynomial currently representing the !! solution, evaluated at `t = tcur`. !! !! * `rwork(lenrw-neq+1:) = acor` = `lenrw-neq+1` array !! of size `neq` used for the accumulated !! corrections on each step, scaled in the output !! to represent the estimated local error in `y` !! on the last step. this is the vector `e` in !! the description of the error control. it is !! defined only on a successful return from [[dvode]]. integer,intent(in) :: neq !! the size of the ode system (number of first order !! ordinary differential equations). !! `neq` may not be increased during the problem, but !! can be decreased (with `istate = 3` in the input). integer,intent(in) :: itol !! an indicator for the type of error control. !! 1 or 2 according as `atol` is a scalar or array. !! see description under `atol`. integer,intent(in) :: itask !! an index specifying the task to be performed. !! input only. itask has the following values and meanings: !! !! * 1 means normal computation of output values of y(t) at !! t = tout (by overshooting and interpolating). !! * 2 means take one step only and return. !! * 3 means stop at the first internal mesh point at or !! beyond t = tout and return. !! * 4 means normal computation of output values of y(t) at !! t = tout but without overshooting t = tcrit. !! tcrit must be input as rwork(1). tcrit may be equal to !! or beyond tout, but not behind it in the direction of !! integration. this option is useful if the problem !! has a singularity at or beyond t = tcrit. !! * 5 means take one step, without passing tcrit, and return. !! tcrit must be input as rwork(1). !! !! note: if itask = 4 or 5 and the solver reaches tcrit !! (within roundoff), it will return t = tcrit (exactly) to !! indicate this (unless itask = 4 and tout comes before tcrit, !! in which case answers at t = tout are returned first). integer,intent(inout) :: istate !! an index used for input and output to specify the !! the state of the calculation. !! !! in the input, the values of istate are as follows: !! !! 1. means this is the first call for the problem !! (initializations will be done). see note below. !! 2. means this is not the first call, and the calculation !! is to continue normally, with no change in any input !! parameters except possibly tout and itask. !! (if itol, rtol, and/or atol are changed between calls !! with istate = 2, the new values will be used but not !! tested for legality.) !! 3. means this is not the first call, and the !! calculation is to continue normally, but with !! a change in input parameters other than !! tout and itask. changes are allowed in !! neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu, !! and any of the optional input except h0. !! (see iwork description for ml and mu.) !! !! note: a preliminary call with tout = t is not counted !! as a first call here, as no initialization or checking of !! input is done. (such a call is sometimes useful to include !! the initial conditions in the output.) !! thus the first call for which tout /= t requires !! istate = 1 in the input. !! !! in the output, istate has the following values and meanings: !! !! * 1 means nothing was done, as tout was equal to t with !! istate = 1 in the input. !! * 2 means the integration was performed successfully. !! * -1 means an excessive amount of work (more than mxstep !! steps) was done on this call, before completing the !! requested task, but the integration was otherwise !! successful as far as t. (mxstep is an optional input !! and is normally 500.) to continue, the user may !! simply reset istate to a value > 1 and call again. !! (the excess work step counter will be reset to 0.) !! in addition, the user may increase mxstep to avoid !! this error return. (see optional input below.) !! * -2 means too much accuracy was requested for the precision !! of the machine being used. this was detected before !! completing the requested task, but the integration !! was successful as far as t. to continue, the tolerance !! parameters must be reset, and istate must be set !! to 3. the optional output tolsf may be used for this !! purpose. (note: if this condition is detected before !! taking any steps, then an illegal input return !! (istate = -3) occurs instead.) !! * -3 means illegal input was detected, before taking any !! integration steps. see written message for details. !! note: if the solver detects an infinite loop of calls !! to the solver with illegal input, it will cause !! the run to stop. !! * -4 means there were repeated error test failures on !! one attempted step, before completing the requested !! task, but the integration was successful as far as t. !! the problem may have a singularity, or the input !! may be inappropriate. !! * -5 means there were repeated convergence test failures on !! one attempted step, before completing the requested !! task, but the integration was successful as far as t. !! this may be caused by an inaccurate jacobian matrix, !! if one is being used. !! * -6 means ewt(i) became zero for some i during the !! integration. pure relative error control (atol(i)=0.0) !! was requested on a variable which has now vanished. !! the integration was successful as far as t. !! !! note: since the normal output value of istate is 2, !! it does not need to be reset for normal continuation. !! also, since a negative input value of istate will be !! regarded as illegal, a negative output value requires the !! user to change it, and possibly other input, before !! calling the solver again. integer,intent(in) :: iopt !! an integer flag to specify whether or not any optional !! input is being used on this call. !! the optional input is listed separately below. !! !! * iopt = 0 means no optional input is being used. !! default values will be used in all cases. !! * iopt = 1 means optional input is being used. integer,intent(in) :: liw !! the length of the array iwork, as declared by the user. !! (this will be checked by the solver.) integer :: iwork(liw) !! an integer work array. the length of iwork must be at least: !! !! * `30` if miter = 0 or 3 (mf = 10, 13, 20, 23), or !! * `30 + neq` otherwise (abs(mf) = 11,12,14,15,21,22,24,25). !! !! the first 30 words of `iwork` are reserved for conditional and !! optional input and optional output. !! !! the following 2 words in `iwork` are conditional input: !! !! * `iwork(1) = ml` !! * `iwork(2) = mu` !! !! these are the lower and upper !! half-bandwidths, respectively, of the !! banded jacobian, excluding the main diagonal. !! the band is defined by the matrix locations !! `(i,j)` with `i-ml <= j <= i+mu`. `ml` and `mu` !! must satisfy `0 <= ml,mu <= neq-1`. !! these are required if `miter` is 4 or 5, and !! ignored otherwise. `ml` and `mu` may in fact be !! the band parameters for a matrix to which !! `df/dy` is only approximately equal. !! !! The following optional input requires `iopt = 1`, and in that !! case all of this input is examined. a value of zero for any !! of these optional input variables will cause the default value to be !! used. thus to use a subset of the optional input, simply preload !! locations 5 to 10 in `iwork` to 0, and !! then set those of interest to nonzero values. !! !! * `iwork(5) = maxord` = the maximum order to be allowed. the default !! value is 12 if `meth = 1`, and 5 if `meth = 2`. !! if `maxord` exceeds the default value, it will !! be reduced to the default value. !! if `maxord` is changed during the problem, it may !! cause the current order to be reduced. !! !! * `iwork(6) = mxstep` = maximum number of (internally defined) steps !! allowed during one call to the solver. !! the default value is 500. !! !! * `iwork(7) = mxhnil` = maximum number of messages printed (per problem) !! warning that `t + h = t` on a step (`h` = step size). !! this must be positive to result in a non-default !! value. the default value is 10. !! !! as optional additional output from dvode, the variables listed !! below are quantities related to the performance of dvode !! which are available to the user. !! except where stated otherwise, all of this output is defined !! on any successful return from dvode, and on any return with !! istate = -1, -2, -4, -5, or -6. on an illegal input return !! (istate = -3), they will be unchanged from their existing values !! (if any), except possibly for `lenrw`, and `leniw`. !! on any error return, output relevant to the error will be defined, !! as noted below: !! !! * `iwork(11) = nst` = the number of steps taken for the problem so far. !! * `iwork(12) = nfe` = the number of f evaluations for the problem so far. !! * `iwork(13) = nje` = the number of jacobian evaluations so far. !! * `iwork(14) = nqu` = the method order last used (successfully). !! * `iwork(15) = nqcur` = the order to be attempted on the next step. !! * `iwork(16) = imxer` = the index of the component of largest magnitude in !! the weighted local error vector ( e(i)/ewt(i) ), !! on an error return with istate = -4 or -5. !! * `iwork(17) = lenrw` = the length of rwork actually required. !! this is defined on normal returns and on an illegal !! input return for insufficient storage. !! * `iwork(18) = leniw` = the length of iwork actually required. !! this is defined on normal returns and on an illegal !! input return for insufficient storage. !! * `iwork(19) = nlu` = the number of matrix lu decompositions so far. !! * `iwork(20) = nni` = the number of nonlinear (newton) iterations so far. !! * `iwork(21) = ncfn` = the number of convergence failures of the nonlinear !! solver so far. !! * `iwork(22) = netf` = the number of error test failures of the integrator !! so far. integer,intent(in) :: mf !! method flag. standard values are: !! !! * 10 for nonstiff (adams) method, no jacobian used. !! * 21 for stiff (bdf) method, user-supplied full jacobian. !! * 22 for stiff method, internally generated full jacobian. !! * 24 for stiff method, user-supplied banded jacobian. !! * 25 for stiff method, internally generated banded jacobian. !! !! the complete set of legal values of !! mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, !! -11, -12, -14, -15, -21, -22, -24, -25. !! !! `mf` is a signed two-digit integer, `mf = jsv*(10*meth + miter)` !! !! `jsv = sign(mf)` indicates the jacobian-saving strategy: !! !! * jsv = 1 means a copy of the jacobian is saved for reuse !! in the corrector iteration algorithm. !! * jsv = -1 means a copy of the jacobian is not saved !! (valid only for miter = 1, 2, 4, or 5). !! !! `meth` indicates the basic linear multistep method: !! !! * meth = 1 means the implicit adams method. !! * meth = 2 means the method based on backward !! differentiation formulas (bdf-s). !! !! `miter` indicates the corrector iteration method: !! !! * miter = 0 means functional iteration (no jacobian matrix !! is involved). !! * miter = 1 means chord iteration with a user-supplied !! full (neq by neq) jacobian. !! * miter = 2 means chord iteration with an internally !! generated (difference quotient) full jacobian !! (using neq extra calls to f per df/dy value). !! * miter = 3 means chord iteration with an internally !! generated diagonal jacobian approximation !! (using 1 extra call to f per df/dy evaluation). !! * miter = 4 means chord iteration with a user-supplied !! banded jacobian. !! * miter = 5 means chord iteration with an internally !! generated banded jacobian (using ml+mu+1 extra !! calls to f per df/dy evaluation). !! !! if miter = 1 or 4, the user must supply a subroutine `jac` !! (the name is arbitrary) as described above under `jac`. !! for other values of miter, a dummy argument can be used. logical :: ihit real(wp) :: atoli , big , ewti , h0 , hmax , hmx , & rh , rtoli , size , tcrit , & tnext , tolsf , tp integer :: i , ier , iflag , imxer , jco , kgo , leniw , lenj , & lenp , lenrw , lenwm , lf0 , mband , mfa , ml , & mu , niter , nslast character(len=80) :: msg integer,parameter :: mxhnl0 = 10 integer,parameter :: mxstp0 = 500 real(wp),parameter :: one = 1.0_wp real(wp),parameter :: two = 2.0_wp real(wp),parameter :: four = 4.0_wp real(wp),parameter :: pt2 = 0.2_wp real(wp),parameter :: hun = 100.0_wp !----------------------------------------------------------------------- ! block a. ! this code block is executed on every call. ! it tests istate and itask for legality and branches appropriately. ! if istate > 1 but the flag init shows that initialization has ! not yet been done, an error return occurs. ! if istate = 1 and tout = t, return immediately. !----------------------------------------------------------------------- if ( istate<1 .or. istate>3 ) then !----------------------------------------------------------------------- ! block i. ! the following block handles all error returns due to illegal input ! (istate = -3), as detected before calling the core integrator. ! first the error message routine is called. if the illegal input ! is a negative istate, the run is aborted (apparent infinite loop). !----------------------------------------------------------------------- msg = 'dvode-- istate (=i1) illegal ' call me%xerrwd(msg,30,1,1,1,istate,0,0,zero,zero) if ( istate>=0 ) then istate = -3 return end if msg = 'dvode-- run aborted: apparent infinite loop ' call me%xerrwd(msg,50,303,2,0,0,0,0,zero,zero) return end if if ( itask<1 .or. itask>5 ) then msg = 'dvode-- itask (=i1) illegal ' call me%xerrwd(msg,30,2,1,1,itask,0,0,zero,zero) istate = -3 return end if if ( istate==1 ) then me%dat%init = 0 if ( tout==t ) return else if ( me%dat%init/=1 ) then msg = 'dvode-- istate (=i1) > 1 but dvode not initialized ' call me%xerrwd(msg,60,3,1,1,istate,0,0,zero,zero) istate = -3 return else if ( istate==2 ) then goto 50 endif !----------------------------------------------------------------------- ! block b. ! the next code block is executed for the initial call (istate = 1), ! or for a continuation call with parameter changes (istate = 3). ! it contains checking of all input and various initializations. ! ! first check legality of the non-optional input neq, itol, iopt, ! mf, ml, and mu. !----------------------------------------------------------------------- if ( neq<=0 ) then msg = 'dvode-- neq (=i1) < 1 ' call me%xerrwd(msg,30,4,1,1,neq,0,0,zero,zero) istate = -3 return end if if ( istate/=1 ) then if ( neq>me%dat%n ) then msg = 'dvode-- istate = 3 and neq increased (i1 to i2) ' call me%xerrwd(msg,50,5,1,2,me%dat%n,neq,0,zero,zero) istate = -3 return endif endif me%dat%n = neq if ( itol<1 .or. itol>4 ) then msg = 'dvode-- itol (=i1) illegal ' call me%xerrwd(msg,30,6,1,1,itol,0,0,zero,zero) istate = -3 return else if ( iopt<0 .or. iopt>1 ) then msg = 'dvode-- iopt (=i1) illegal ' call me%xerrwd(msg,30,7,1,1,iopt,0,0,zero,zero) istate = -3 return else me%dat%jsv = sign(1,mf) mfa = abs(mf) me%dat%meth = mfa/10 me%dat%miter = mfa - 10*me%dat%meth if (( me%dat%meth<1 .or. me%dat%meth>2 ) .or. & ( me%dat%miter<0 .or. me%dat%miter>5 )) then msg = 'dvode-- mf (=i1) illegal ' call me%xerrwd(msg,30,8,1,1,mf,0,0,zero,zero) istate = -3 return end if if ( me%dat%miter>3 ) then ml = iwork(1) mu = iwork(2) if ( ml<0 .or. ml>=me%dat%n ) then msg = 'dvode-- ml (=i1) illegal: <0 or >=neq (=i2)' call me%xerrwd(msg,50,9,1,2,ml,neq,0,zero,zero) istate = -3 return else if ( mu<0 .or. mu>=me%dat%n ) then msg = 'dvode-- mu (=i1) illegal: <0 or >=neq (=i2)' call me%xerrwd(msg,50,10,1,2,mu,neq,0,zero,zero) istate = -3 return endif endif ! next process and check the optional input. --------------------------- if ( iopt==1 ) then me%dat%maxord = iwork(5) if ( me%dat%maxord<0 ) then msg = 'dvode-- maxord (=i1) < 0 ' call me%xerrwd(msg,30,11,1,1,me%dat%maxord,0,0,zero,zero) istate = -3 return else if ( me%dat%maxord==0 ) me%dat%maxord = 100 me%dat%maxord = min(me%dat%maxord,mord(me%dat%meth)) me%dat%mxstep = iwork(6) if ( me%dat%mxstep<0 ) then msg = 'dvode-- mxstep (=i1) < 0 ' call me%xerrwd(msg,30,12,1,1,me%dat%mxstep,0,0,zero,zero) istate = -3 return else if ( me%dat%mxstep==0 ) me%dat%mxstep = mxstp0 me%dat%mxhnil = iwork(7) if ( me%dat%mxhnil<0 ) then msg = 'dvode-- mxhnil (=i1) < 0 ' call me%xerrwd(msg,30,13,1,1,me%dat%mxhnil,0,0,zero,zero) istate = -3 return else if ( me%dat%mxhnil==0 ) me%dat%mxhnil = mxhnl0 if ( istate==1 ) then h0 = rwork(5) if ( (tout-t)*h0<zero ) then msg = 'dvode-- tout (=r1) behind t (=r2) ' call me%xerrwd(msg,40,14,1,0,0,0,2,tout,t) msg = ' integration direction is given by h0 (=r1) ' call me%xerrwd(msg,50,14,1,0,0,0,1,h0,zero) istate = -3 return endif endif hmax = rwork(6) if ( hmax<zero ) then msg = 'dvode-- hmax (=r1) < 0.0 ' call me%xerrwd(msg,30,15,1,0,0,0,1,hmax,zero) istate = -3 return else me%dat%hmxi = zero if ( hmax>zero ) me%dat%hmxi = one/hmax me%dat%hmin = rwork(7) if ( me%dat%hmin<zero ) then msg = 'dvode-- hmin (=r1) < 0.0 ' call me%xerrwd(msg,30,16,1,0,0,0,1,me%dat%hmin,zero) istate = -3 return endif endif endif endif endif else me%dat%maxord = mord(me%dat%meth) me%dat%mxstep = mxstp0 me%dat%mxhnil = mxhnl0 if ( istate==1 ) h0 = zero me%dat%hmxi = zero me%dat%hmin = zero endif !----------------------------------------------------------------------- ! set work array pointers and check lengths lrw and liw. ! pointers to segments of rwork and iwork are named by prefixing l to ! the name of the segment. e.g., the segment yh starts at rwork(lyh). ! segments of rwork (in order) are denoted yh, wm, ewt, savf, acor. ! within wm, locjs is the location of the saved jacobian (jsv > 0). !----------------------------------------------------------------------- me%dat%lyh = 21 if ( istate==1 ) me%dat%nyh = me%dat%n me%dat%lwm = me%dat%lyh + (me%dat%maxord+1)*me%dat%nyh jco = max(0,me%dat%jsv) select case (me%dat%miter) case(0) lenwm = 0 case(1:2) lenwm = 2 + (1+jco)*me%dat%n*me%dat%n me%dat%locjs = me%dat%n*me%dat%n + 3 case(3) lenwm = 2 + me%dat%n case(4:5) mband = ml + mu + 1 lenp = (mband+ml)*me%dat%n lenj = mband*me%dat%n lenwm = 2 + lenp + jco*lenj me%dat%locjs = lenp + 3 end select me%dat%lewt = me%dat%lwm + lenwm me%dat%lsavf = me%dat%lewt + me%dat%n me%dat%lacor = me%dat%lsavf + me%dat%n lenrw = me%dat%lacor + me%dat%n - 1 iwork(17) = lenrw me%dat%liwm = 1 leniw = 30 + me%dat%n if ( me%dat%miter==0 .or. me%dat%miter==3 ) leniw = 30 iwork(18) = leniw if ( lenrw>lrw ) then msg = 'dvode-- rwork length needed, lenrw (=i1), exceeds lrw (=i2)' call me%xerrwd(msg,60,17,1,2,lenrw,lrw,0,zero,zero) istate = -3 return else if ( leniw>liw ) then msg = 'dvode-- iwork length needed, leniw (=i1), exceeds liw (=i2)' call me%xerrwd(msg,60,18,1,2,leniw,liw,0,zero,zero) istate = -3 return else ! check rtol and atol for legality. ------------------------------------ rtoli = rtol(1) atoli = atol(1) do i = 1 , me%dat%n if ( itol>=3 ) rtoli = rtol(i) if ( itol==2 .or. itol==4 ) atoli = atol(i) if ( rtoli<zero ) then msg = 'dvode-- rtol(i1) is r1 < 0.0 ' call me%xerrwd(msg,40,19,1,1,i,0,1,rtoli,zero) istate = -3 return end if if ( atoli<zero ) then msg = 'dvode-- atol(i1) is r1 < 0.0 ' call me%xerrwd(msg,40,20,1,1,i,0,1,atoli,zero) istate = -3 return end if enddo if ( istate==1 ) then !----------------------------------------------------------------------- ! block c. ! the next block is for the initial call only (istate = 1). ! it contains all remaining initializations, the initial call to f, ! and the calculation of the initial step size. ! the error weights in ewt are inverted after being loaded. !----------------------------------------------------------------------- me%dat%uround = epmach me%dat%tn = t if ( itask==4 .or. itask==5 ) then tcrit = rwork(1) if ( (tcrit-tout)*(tout-t)<zero ) then msg = 'dvode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ' call me%xerrwd(msg,60,25,1,0,0,0,2,tcrit,tout) istate = -3 return end if if ( h0/=zero .and. (t+h0-tcrit)*h0>zero ) & h0 = tcrit - t endif me%dat%jstart = 0 if ( me%dat%miter>0 ) rwork(me%dat%lwm) = sqrt(me%dat%uround) me%dat%ccmxj = pt2 me%dat%msbj = 50 me%dat%nhnil = 0 me%dat%nst = 0 me%dat%nje = 0 me%dat%nni = 0 me%dat%ncfn = 0 me%dat%netf = 0 me%dat%nlu = 0 me%dat%nslj = 0 nslast = 0 me%dat%hu = zero me%dat%nqu = 0 ! initial call to f. (lf0 points to yh(*,2).) ------------------------- lf0 = me%dat%lyh + me%dat%nyh call me%f(me%dat%n,t,y(1:me%dat%n),rwork(lf0)) me%dat%nfe = 1 ! load the initial value vector in yh. --------------------------------- call dcopy(me%dat%n,y,1,rwork(me%dat%lyh),1) ! load and invert the ewt array. (h is temporarily set to 1.0.) ------- me%dat%nq = 1 me%dat%h = one call me%dewset(me%dat%n,itol,rtol(1:me%dat%n),atol(1:me%dat%n),& rwork(me%dat%lyh), rwork(me%dat%lewt)) do i = 1 , me%dat%n if ( rwork(i+me%dat%lewt-1)<=zero ) then ewti = rwork(me%dat%lewt+i-1) msg = 'dvode-- ewt(i1) is r1 <= 0.0 ' call me%xerrwd(msg,40,21,1,1,i,0,1,ewti,zero) istate = -3 return end if rwork(i+me%dat%lewt-1) = one/rwork(i+me%dat%lewt-1) enddo if ( h0==zero ) then ! call dvhin to set initial step size h0 to be attempted. -------------- call me%dvhin(me%dat%n,t,rwork(me%dat%lyh),rwork(lf0), & tout,me%dat%uround,rwork(me%dat%lewt),itol,& atol,y,rwork(me%dat%lacor),h0,niter,ier) me%dat%nfe = me%dat%nfe + niter if ( ier/=0 ) then msg = & 'dvode-- tout (=r1) too close to t(=r2) to start integration' call me%xerrwd(msg,60,22,1,0,0,0,2,tout,t) istate = -3 return endif endif ! adjust h0 if necessary to meet hmax bound. --------------------------- rh = abs(h0)*me%dat%hmxi if ( rh>one ) h0 = h0/rh ! load h with h0 and scale yh(*,2) by h0. ------------------------------ me%dat%h = h0 call dscal(me%dat%n,h0,rwork(lf0),1) goto 200 else ! if istate = 3, set flag to signal parameter changes to dvstep. ------- me%dat%jstart = -1 ! maxord was reduced below nq. copy yh(*,maxord+2) into savf. --------- if ( me%dat%nq>me%dat%maxord ) & call dcopy(me%dat%n,rwork(me%dat%lwm),1,rwork(me%dat%lsavf),1) ! reload wm(1) = rwork(lwm), since lwm may have changed. --------------- if ( me%dat%miter>0 ) rwork(me%dat%lwm) = sqrt(me%dat%uround) endif endif endif !----------------------------------------------------------------------- ! block d. ! the next code block is for continuation calls only (istate = 2 or 3) ! and is to check stop conditions before taking a step. !----------------------------------------------------------------------- 50 nslast = me%dat%nst me%dat%kuth = 0 select case (itask) case (2) goto 100 case (3) tp = me%dat%tn - me%dat%hu*(one+hun*me%dat%uround) if ( (tp-tout)*me%dat%h>zero ) then msg = 'dvode-- itask = i1 and tout (=r1) behind tcur - hu (= r2) ' call me%xerrwd(msg,60,23,1,1,itask,0,2,tout,tp) istate = -3 return else if ( (me%dat%tn-tout)*me%dat%h>=zero ) then call set_y_and_t() if ( itask==4 .or. itask==5 ) then if ( ihit ) t = tcrit endif call continue() return end if goto 100 endif case (4) tcrit = rwork(1) if ( (me%dat%tn-tcrit)*me%dat%h>zero ) then msg = 'dvode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ' call me%xerrwd(msg,60,24,1,0,0,0,2,tcrit,me%dat%tn) istate = -3 return end if if ( (tcrit-tout)*me%dat%h<zero ) then msg = 'dvode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ' call me%xerrwd(msg,60,25,1,0,0,0,2,tcrit,tout) istate = -3 return end if if ( (me%dat%tn-tout)*me%dat%h>=zero ) then call me%dvindy(tout,0,rwork(me%dat%lyh),me%dat%nyh,y,iflag) if ( iflag/=0 ) then msg = 'dvode-- trouble from dvindy. itask = i1, tout = r1. ' call me%xerrwd(msg,60,27,1,1,itask,0,1,tout,zero) istate = -3 return end if t = tout call continue() return endif case (5) tcrit = rwork(1) if ( (me%dat%tn-tcrit)*me%dat%h>zero ) then msg = 'dvode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ' call me%xerrwd(msg,60,24,1,0,0,0,2,tcrit,me%dat%tn) istate = -3 return end if case default if ( (me%dat%tn-tout)*me%dat%h<zero ) goto 100 call me%dvindy(tout,0,rwork(me%dat%lyh),me%dat%nyh,y,iflag) if ( iflag/=0 ) then msg = 'dvode-- trouble from dvindy. itask = i1, tout = r1. ' call me%xerrwd(msg,60,27,1,1,itask,0,1,tout,zero) istate = -3 return end if t = tout call continue() return end select hmx = abs(me%dat%tn) + abs(me%dat%h) ihit = abs(me%dat%tn-tcrit)<=hun*me%dat%uround*hmx if ( ihit ) then call set_y_and_t() if ( itask==4 .or. itask==5 ) then if ( ihit ) t = tcrit endif call continue() return end if tnext = me%dat%tn + me%dat%hnew*(one+four*me%dat%uround) if ( (tnext-tcrit)*me%dat%h>zero ) then me%dat%h = (tcrit-me%dat%tn)*(one-four*me%dat%uround) me%dat%kuth = 1 endif !----------------------------------------------------------------------- ! block e. ! the next block is normally executed for all calls and contains ! the call to the one-step core integrator dvstep. ! ! this is a looping point for the integration steps. ! ! first check for too many steps being taken, update ewt (if not at ! start of problem), check for too much accuracy being requested, and ! check for h below the roundoff level in t. !----------------------------------------------------------------------- 100 if ( (me%dat%nst-nslast)>=me%dat%mxstep ) then !----------------------------------------------------------------------- ! block h. ! the following block handles all unsuccessful returns other than ! those for illegal input. first the error message routine is called. ! if there was an error test or convergence test failure, imxer is set. ! then y is loaded from yh, and t is set to tn. ! the optional output is loaded into the work arrays before returning. !----------------------------------------------------------------------- ! the maximum number of steps was taken before reaching tout. ---------- msg = 'dvode-- at current t (=r1), mxstep (=i1) steps ' call me%xerrwd(msg,50,201,1,0,0,0,0,zero,zero) msg = ' taken on this call before reaching tout ' call me%xerrwd(msg,50,201,1,1,me%dat%mxstep,0,1,me%dat%tn,zero) istate = -1 call finish() return else call me%dewset(me%dat%n,itol,rtol(1:me%dat%n),atol(1:me%dat%n),& rwork(me%dat%lyh),rwork(me%dat%lewt)) do i = 1 , me%dat%n if ( rwork(i+me%dat%lewt-1)<=zero ) then call ewt_error() call finish() return end if rwork(i+me%dat%lewt-1) = one/rwork(i+me%dat%lewt-1) enddo endif 200 tolsf = me%dat%uround*me%dvnorm(me%dat%n,rwork(me%dat%lyh),rwork(me%dat%lewt)) if ( tolsf<=one ) then if ( (me%dat%tn+me%dat%h)==me%dat%tn ) then me%dat%nhnil = me%dat%nhnil + 1 if ( me%dat%nhnil<=me%dat%mxhnil ) then msg = 'dvode-- warning: internal t (=r1) and h (=r2) are' call me%xerrwd(msg,50,101,1,0,0,0,0,zero,zero) msg = ' such that in the machine, t + h = t on the next step ' call me%xerrwd(msg,60,101,1,0,0,0,0,zero,zero) msg = ' (h = step size). solver will continue anyway' call me%xerrwd(msg,50,101,1,0,0,0,2,me%dat%tn,me%dat%h) if ( me%dat%nhnil>=me%dat%mxhnil ) then msg = 'dvode-- above warning has been issued i1 times. ' call me%xerrwd(msg,50,102,1,0,0,0,0,zero,zero) msg = ' it will not be issued again for this problem' call me%xerrwd(msg,50,102,1,1,me%dat%mxhnil,0,0,zero,zero) endif endif endif call me%dvstep(y,& ! y rwork(me%dat%lyh),& ! yh me%dat%nyh,& ! nyh rwork(me%dat%lyh),& ! yh rwork(me%dat%lewt), & ! ewt rwork(me%dat%lsavf),& ! savf y,& ! vsav rwork(me%dat%lacor),& ! acor rwork(me%dat%lwm),& ! wm iwork(me%dat%liwm)) ! iwm kgo = 1 - me%dat%kflag ! branch on kflag. note: in this version, kflag can not be set to -3. ! kflag == 0, -1, -2 select case (kgo) case (2) ! kflag = -1. error test failed repeatedly or with abs(h) = hmin. ----- msg = 'dvode-- at t(=r1) and step size h(=r2), the error' call me%xerrwd(msg,50,204,1,0,0,0,0,zero,zero) msg = ' test failed repeatedly or with abs(h) = hmin' call me%xerrwd(msg,50,204,1,0,0,0,2,me%dat%tn,me%dat%h) istate = -4 call compute_imxer() call finish() return case (3) ! kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ---- msg = 'dvode-- at t (=r1) and step size h (=r2), the ' call me%xerrwd(msg,50,205,1,0,0,0,0,zero,zero) msg = ' corrector convergence failed repeatedly ' call me%xerrwd(msg,50,205,1,0,0,0,0,zero,zero) msg = ' or with abs(h) = hmin ' call me%xerrwd(msg,30,205,1,0,0,0,2,me%dat%tn,me%dat%h) istate = -5 call compute_imxer() call finish() return case default !----------------------------------------------------------------------- ! block f. ! the following block handles the case of a successful return from the ! core integrator (kflag = 0). test for stop conditions. !----------------------------------------------------------------------- me%dat%init = 1 me%dat%kuth = 0 select case (itask) case (2) case (3) ! itask = 3. jump to exit if tout was reached. ------------------------ if ( (me%dat%tn-tout)*me%dat%h<zero ) goto 100 case (4) ! itask = 4. see if tout or tcrit was reached. adjust h if necessary. if ( (me%dat%tn-tout)*me%dat%h<zero ) then hmx = abs(me%dat%tn) + abs(me%dat%h) ihit = abs(me%dat%tn-tcrit)<=hun*me%dat%uround*hmx if ( .not.(ihit) ) then tnext = me%dat%tn + me%dat%hnew*(one+four*me%dat%uround) if ( (tnext-tcrit)*me%dat%h>zero ) then me%dat%h = (tcrit-me%dat%tn)*(one-four*me%dat%uround) me%dat%kuth = 1 endif goto 100 endif else call me%dvindy(tout,0,rwork(me%dat%lyh),me%dat%nyh,y,iflag) t = tout call continue() return endif case (5) ! itask = 5. see if tcrit was reached and jump to exit. --------------- hmx = abs(me%dat%tn) + abs(me%dat%h) ihit = abs(me%dat%tn-tcrit)<=hun*me%dat%uround*hmx case default ! itask = 1. if tout has been reached, interpolate. ------------------- if ( (me%dat%tn-tout)*me%dat%h<zero ) goto 100 call me%dvindy(tout,0,rwork(me%dat%lyh),me%dat%nyh,y,iflag) t = tout call continue() return end select end select else tolsf = tolsf*two if ( me%dat%nst==0 ) then msg = 'dvode-- at start of problem, too much accuracy ' call me%xerrwd(msg,50,26,1,0,0,0,0,zero,zero) msg = ' requested for precision of machine: see tolsf (=r1) ' call me%xerrwd(msg,60,26,1,0,0,0,1,tolsf,zero) rwork(14) = tolsf istate = -3 return else ! too much accuracy requested for machine precision. ------------------- msg = 'dvode-- at t (=r1), too much accuracy requested ' call me%xerrwd(msg,50,203,1,0,0,0,0,zero,zero) msg = ' for precision of machine: see tolsf (=r2) ' call me%xerrwd(msg,50,203,1,0,0,0,2,me%dat%tn,tolsf) rwork(14) = tolsf istate = -2 call finish() return endif endif call set_y_and_t() if ( itask==4 .or. itask==5 ) then if ( ihit ) t = tcrit endif call continue() contains !----------------------------------------------------------------------- ! block g. ! the following block handles all successful returns from dvode. ! if itask /= 1, y is loaded from yh and t is set accordingly. ! istate is set to 2, and the optional output is loaded into the work ! arrays before returning. !----------------------------------------------------------------------- subroutine set_y_and_t() !! set y vector, t call dcopy(me%dat%n,rwork(me%dat%lyh),1,y,1) t = me%dat%tn end subroutine set_y_and_t subroutine continue() istate = 2 rwork(11) = me%dat%hu rwork(12) = me%dat%hnew rwork(13) = me%dat%tn iwork(11) = me%dat%nst iwork(12) = me%dat%nfe iwork(13) = me%dat%nje iwork(14) = me%dat%nqu iwork(15) = me%dat%newq iwork(19) = me%dat%nlu iwork(20) = me%dat%nni iwork(21) = me%dat%ncfn iwork(22) = me%dat%netf end subroutine continue subroutine ewt_error() !! ewt(i) <= 0.0 for some i (not at start of problem). ewti = rwork(me%dat%lewt+i-1) msg = 'dvode-- at t (=r1), ewt(i1) has become r2 <= 0.' call me%xerrwd(msg,50,202,1,1,i,0,2,me%dat%tn,ewti) istate = -6 end subroutine ewt_error subroutine compute_imxer() !! compute imxer if relevant. big = zero imxer = 1 do i = 1 , me%dat%n size = abs(rwork(i+me%dat%lacor-1)*rwork(i+me%dat%lewt-1)) if ( big<size ) then big = size imxer = i endif enddo iwork(16) = imxer end subroutine compute_imxer subroutine finish() !! set y vector, t, and optional output. call set_y_and_t() rwork(11) = me%dat%hu rwork(12) = me%dat%h rwork(13) = me%dat%tn iwork(11) = me%dat%nst iwork(12) = me%dat%nfe iwork(13) = me%dat%nje iwork(14) = me%dat%nqu iwork(15) = me%dat%nq iwork(19) = me%dat%nlu iwork(20) = me%dat%nni iwork(21) = me%dat%ncfn iwork(22) = me%dat%netf end subroutine finish end subroutine dvode