dstodpk Subroutine

public subroutine dstodpk(Neq, Y, Yh, Nyh, Yh1, Ewt, Savf, Savx, Acor, Wm, Iwm, f, jac, psol)

DSTODPK performs one step of the integration of an initial value problem for a system of Ordinary Differential Equations.

The following changes were made to generate Subroutine DSTODPK from Subroutine DSTODE:

1. The array SAVX was added to the call sequence.
2. PJAC and SLVS were replaced by PSOL in the call sequence.
3. The Common block /DLPK01/ was added for communication.
4. The test constant EPCON is loaded into Common below statement numbers 125 and 155, and used below statement 400.
5. The Newton iteration counter MNEWT is set below 220 and 400.
6. The call to PJAC was replaced with a call to DPKSET (fixed name), with a longer call sequence, called depending on JACFLG.
7. The corrector residual is stored in SAVX (not Y) at 360, and the solution vector is in SAVX in the 380 loop.
8. SLVS was renamed DSOLPK and includes NEQ, SAVX, EWT, F, and JAC. SAVX was added because DSOLPK now needs Y and SAVF undisturbed.
9. The nonlinear convergence failure count NCFN is set at 430.

Note: DSTODPK is independent of the value of the iteration method indicator MITER, when this is .ne. 0, and hence is independent of the type of chord method used, or the Jacobian structure.

Communication with DSTODPK is done with the following variables:

NEQ

integer array containing problem size in NEQ(1), and passed as the NEQ argument in all calls to F and JAC.

Y

an array of length .ge. N used as the Y argument in all calls to F and JAC.

YH

an NYH by LMAX array containing the dependent variables and their approximate scaled derivatives, where LMAX = MAXORD + 1. YH(i,j+1) contains the approximate j-th derivative of y(i), scaled by H**j/factorial(j) (j = 0,1,…,NQ). On entry for the first step, the first two columns of YH must be set from the initial values.

NYH

a constant integer .ge. N, the first dimension of YH.

YH1

a one-dimensional array occupying the same space as YH.

EWT

an array of length N containing multiplicative weights for local error measurements. Local errors in y(i) are compared to 1.0/EWT(i) in various error tests.

SAVF

an array of working storage, of length N. Also used for input of YH(*,MAXORD+2) when JSTART = -1 and MAXORD .lt. the current order NQ.

SAVX

an array of working storage, of length N.

ACOR

a work array of length N, used for the accumulated corrections. On a successful return, ACOR(i) contains the estimated one-step local error in y(i).

WM,IWM

real and integer work arrays associated with matrix operations in chord iteration (MITER .ne. 0).

CCMAX

maximum relative change in H*EL0 before DPKSET is called.

H

the step size to be attempted on the next step. H is altered by the error control algorithm during the problem. H can be either positive or negative, but its sign must remain constant throughout the problem.

HMIN

the minimum absolute value of the step size H to be used.

HMXI

inverse of the maximum absolute value of H to be used. HMXI = 0.0 is allowed and corresponds to an infinite HMAX. HMIN and HMXI may be changed at any time, but will not take effect until the next change of H is considered.

TN

the independent variable. TN is updated on each step taken.

JSTART

an integer used for input only, with the following values and meanings:

     0  perform the first step.
 .gt.0  take a new step continuing from the last.
    -1  take the next step with a new value of H, MAXORD,
        N, METH, MITER, and/or matrix parameters.
    -2  take the next step with a new value of H,
        but with other inputs unchanged.

On return, JSTART is set to 1 to facilitate continuation.

KFLAG

a completion code with the following meanings:

     0  the step was succesful.
    -1  the requested error could not be achieved.
    -2  corrector convergence could not be achieved.
    -3  fatal error in DPKSET or DSOLPK.

A return with KFLAG = -1 or -2 means either ABS(H) = HMIN or 10 consecutive failures occurred. On a return with KFLAG negative, the values of TN and the YH array are as of the beginning of the last step, and H is the last step size attempted.

MAXORD

the maximum order of integration method to be allowed.

MAXCOR

the maximum number of corrector iterations allowed.

MSBP

maximum number of steps between DPKSET calls (MITER .gt. 0).

MXNCF

maximum number of convergence failures allowed.

METH/MITER

the method flags. See description in driver.

N

the number of first-order differential equations.

Arguments

Type	Intent	Optional	Attributes		Name
integer,			dimension(*)	::	Neq
real(kind=dp),			dimension(*)	::	Y
real(kind=dp),	intent(inout),		dimension(Nyh,*)	::	Yh
integer,	intent(in)			::	Nyh
real(kind=dp),	intent(inout),		dimension(*)	::	Yh1
real(kind=dp),			dimension(*)	::	Ewt
real(kind=dp),	intent(inout),		dimension(*)	::	Savf
real(kind=dp),	intent(inout),		dimension(*)	::	Savx
real(kind=dp),	intent(inout),		dimension(*)	::	Acor
real(kind=dp),			dimension(*)	::	Wm
integer,			dimension(*)	::	Iwm
real				::	f
integer				::	jac
real				::	psol