| [n,t,y,f,nstate,tout,ntask,nroot,eps,ewt,ierror,mint,miter,impl,ml,mu,mxord,hmax,work,lenw,iwork,leniw,jacobn,fa,nde,mxstep,g,users,ierflg]=sdriv3(n,t,y,f,nstate,tout,ntask,nroot,eps,ewt,ierror,mint,miter,impl,ml,mu,mxord,hmax,work,lenw,iwork,leniw,jacobn |
function [n,t,y,f,nstate,tout,ntask,nroot,eps,ewt,ierror,mint,miter,impl,ml,mu,mxord,hmax,work,lenw,iwork,leniw,jacobn,fa,nde,mxstep,g,users,ierflg]=sdriv3(n,t,y,f,nstate,tout,ntask,nroot,eps,ewt,ierror,mint,miter,impl,ml,mu,mxord,hmax,work,lenw,iwork,leniw,jacobn,fa,nde,mxstep,g,users,ierflg);
%***BEGIN PROLOGUE SDRIV3
%***PURPOSE The function of SDRIV3 is to solve N ordinary differential
% equations of the form dY(I)/dT = F(Y(I),T), given the
% initial conditions Y(I) = YI. The program has options to
% allow the solution of both stiff and non-stiff differential
% equations. Other important options are available. SDRIV3
% uses single precision arithmetic.
%***LIBRARY SLATEC (SDRIVE)
%***CATEGORY I1A2, I1A1B
%***TYPE SINGLE PRECISION (SDRIV3-S, DDRIV3-D, CDRIV3-C)
%***KEYWORDS GEAR'S METHOD, INITIAL VALUE PROBLEMS, ODE,
% ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, SINGLE PRECISION,
% STIFF
%***AUTHOR Kahaner, D. K., (NIST)
% National Institute of Standards and Technology
% Gaithersburg, MD 20899
% Sutherland, C. D., (LANL)
% Mail Stop D466
% Los Alamos National Laboratory
% Los Alamos, NM 87545
%***DESCRIPTION
%
% I. ABSTRACT .......................................................
%
% The primary function of SDRIV3 is to solve N ordinary differential
% equations of the form dY(I)/dT = F(Y(I),T), given the initial
% conditions Y(I) = YI. The program has options to allow the
% solution of both stiff and non-stiff differential equations. In
% addition, SDRIV3 may be used to solve:
% 1. The initial value problem, A*dY(I)/dT = F(Y(I),T), where A is
% a non-singular matrix depending on Y and T.
% 2. The hybrid differential/algebraic initial value problem,
% A*dY(I)/dT = F(Y(I),T), where A is a vector (whose values may
% depend upon Y and T) some of whose components will be zero
% corresponding to those equations which are algebraic rather
% than differential.
% SDRIV3 is to be called once for each output point of T.
%
% II. PARAMETERS ....................................................
%
% The user should use parameter names in the call sequence of SDRIV3
% for those quantities whose value may be altered by SDRIV3. The
% parameters in the call sequence are:
%
% N = (Input) The number of dependent functions whose solution
% is desired. N must not be altered during a problem.
%
% T = The independent variable. On input for the first call, T
% is the initial point. On output, T is the point at which
% the solution is given.
%
% Y = The vector of dependent variables. Y is used as input on
% the first call, to set the initial values. On output, Y
% is the computed solution vector. This array Y is passed
% in the call sequence of the user-provided routines F,
% JACOBN, FA, USERS, and G. Thus parameters required by
% those routines can be stored in this array in components
% N+1 and above. (Note: Changes by the user to the first
% N components of this array will take effect only after a
% restart, i.e., after setting NSTATE to 1 .)
%
% F = A subroutine supplied by the user. The name must be
% declared EXTERNAL in the user's calling program. This
% subroutine is of the form:
% subroutine F (N, T, Y, YDOT)
% REAL Y(*), YDOT(*)
% .
% .
% YDOT(1) = ...
% .
% .
% YDOT(N) = ...
% end (Sample)
% This computes YDOT = F(Y,T), the right hand side of the
% differential equations. Here Y is a vector of length at
% least N. The actual length of Y is determined by the
% user's declaration in the program which calls SDRIV3.
% Thus the dimensioning of Y in F, while required by FORTRAN
% convention, does not actually allocate any storage. When
% this subroutine is called, the first N components of Y are
% intermediate approximations to the solution components.
% The user should not alter these values. Here YDOT is a
% vector of length N. The user should only compute YDOT(I)
% for I from 1 to N. Normally a return from F passes
% control back to SDRIV3. However, if the user would like
% to abort the calculation, i.e., return control to the
% program which calls SDRIV3, he should set N to zero.
% SDRIV3 will signal this by returning a value of NSTATE
% equal to 6 . Altering the value of N in F has no effect
% on the value of N in the call sequence of SDRIV3.
%
% NSTATE = An integer describing the status of integration. The
% meaning of NSTATE is as follows:
% 1 (Input) Means the first call to the routine. This
% value must be set by the user. On all subsequent
% calls the value of NSTATE should be tested by the
% user, but must not be altered. (As a convenience to
% the user who may wish to put out the initial
% conditions, SDRIV3 can be called with NSTATE=1, and
% TOUT=T. In this case the program will return with
% NSTATE unchanged, i.e., NSTATE=1.)
% 2 (Output) Means a successful integration. If a normal
% continuation is desired (i.e., a further integration
% in the same direction), simply advance TOUT and call
% again. All other parameters are automatically set.
% 3 (Output,Unsuccessful) Means the integrator has taken
% MXSTEP steps without reaching TOUT. The user can
% continue the integration by simply calling SDRIV3
% again.
% 4 (Output,Unsuccessful) Means too much accuracy has
% been requested. EPS has been increased to a value
% the program estimates is appropriate. The user can
% continue the integration by simply calling SDRIV3
% again.
% 5 (Output) A root was found at a point less than TOUT.
% The user can continue the integration toward TOUT by
% simply calling SDRIV3 again.
% 6 (Output,Unsuccessful) N has been set to zero in
% subroutine F.
% 7 (Output,Unsuccessful) N has been set to zero in
% function G. See description of G below.
% 8 (Output,Unsuccessful) N has been set to zero in
% subroutine JACOBN. See description of JACOBN below.
% 9 (Output,Unsuccessful) N has been set to zero in
% subroutine FA. See description of FA below.
% 10 (Output,Unsuccessful) N has been set to zero in
% subroutine USERS. See description of USERS below.
% 11 (Output,Successful) For NTASK = 2 or 3, T is beyond
% TOUT. The solution was obtained by interpolation.
% The user can continue the integration by simply
% advancing TOUT and calling SDRIV3 again.
% 12 (Output,Unsuccessful) The solution could not be
% obtained. The value of IERFLG (see description
% below) for a 'Recoverable' situation indicates the
% type of difficulty encountered: either an illegal
% value for a parameter or an inability to continue the
% solution. For this condition the user should take
% corrective action and reset NSTATE to 1 before
% calling SDRIV3 again. Otherwise the program will
% terminate the run.
%
% TOUT = (Input) The point at which the solution is desired. The
% position of TOUT relative to T on the first call
% determines the direction of integration.
%
% NTASK = (Input) An index specifying the manner of returning the
% solution, according to the following:
% NTASK = 1 Means SDRIV3 will integrate past TOUT and
% interpolate the solution. This is the most
% efficient mode.
% NTASK = 2 Means SDRIV3 will return the solution after
% each internal integration step, or at TOUT,
% whichever comes first. In the latter case,
% the program integrates exactly to TOUT.
% NTASK = 3 Means SDRIV3 will adjust its internal step to
% reach TOUT exactly (useful if a singularity
% exists beyond TOUT.)
%
% NROOT = (Input) The number of equations whose roots are desired.
% If NROOT is zero, the root search is not active. This
% option is useful for obtaining output at points which are
% not known in advance, but depend upon the solution, e.g.,
% when some solution component takes on a specified value.
% The root search is carried out using the user-written
% function G (see description of G below.) SDRIV3 attempts
% to find the value of T at which one of the equations
% changes sign. SDRIV3 can find at most one root per
% equation per internal integration step, and will then
% return the solution either at TOUT or at a root, whichever
% occurs first in the direction of integration. The initial
% point is never reported as a root. The index of the
% equation whose root is being reported is stored in the
% sixth element of IWORK.
% NOTE: NROOT is never altered by this program.
%
% EPS = On input, the requested relative accuracy in all solution
% components. EPS = 0 is allowed. On output, the adjusted
% relative accuracy if the input value was too small. The
% value of EPS should be set as large as is reasonable,
% because the amount of work done by SDRIV3 increases as EPS
% decreases.
%
% EWT = (Input) Problem zero, i.e., the smallest, nonzero,
% physically meaningful value for the solution. (Array,
% possibly of length one. See following description of
% IERROR.) Setting EWT smaller than necessary can adversely
% affect the running time.
%
% IERROR = (Input) Error control indicator. A value of 3 is
% suggested for most problems. Other choices and detailed
% explanations of EWT and IERROR are given below for those
% who may need extra flexibility.
%
% These last three input quantities EPS, EWT and IERROR
% control the accuracy of the computed solution. EWT and
% IERROR are used internally to compute an array YWT. One
% step error estimates divided by YWT(I) are kept less than
% EPS in root mean square norm.
% IERROR (Set by the user) =
% 1 Means YWT(I) = 1. (Absolute error control)
% EWT is ignored.
% 2 Means YWT(I) = ABS(Y(I)), (Relative error control)
% EWT is ignored.
% 3 Means YWT(I) = MAX(ABS(Y(I)), EWT(1)).
% 4 Means YWT(I) = MAX(ABS(Y(I)), EWT(I)).
% This choice is useful when the solution components
% have differing scales.
% 5 Means YWT(I) = EWT(I).
% If IERROR is 3, EWT need only be dimensioned one.
% If IERROR is 4 or 5, the user must dimension EWT at least
% N, and set its values.
%
% MINT = (Input) The integration method indicator.
% MINT = 1 Means the Adams methods, and is used for
% non-stiff problems.
% MINT = 2 Means the stiff methods of Gear (i.e., the
% backward differentiation formulas), and is
% used for stiff problems.
% MINT = 3 Means the program dynamically selects the
% Adams methods when the problem is non-stiff
% and the Gear methods when the problem is
% stiff. When using the Adams methods, the
% program uses a value of MITER=0; when using
% the Gear methods, the program uses the value
% of MITER provided by the user. Only a value
% of IMPL = 0 and a value of MITER = 1, 2, 4, or
% 5 is allowed for this option. The user may
% not alter the value of MINT or MITER without
% restarting, i.e., setting NSTATE to 1.
%
% MITER = (Input) The iteration method indicator.
% MITER = 0 Means functional iteration. This value is
% suggested for non-stiff problems.
% MITER = 1 Means chord method with analytic Jacobian.
% In this case, the user supplies subroutine
% JACOBN (see description below).
% MITER = 2 Means chord method with Jacobian calculated
% internally by finite differences.
% MITER = 3 Means chord method with corrections computed
% by the user-written routine USERS (see
% description of USERS below.) This option
% allows all matrix algebra and storage
% decisions to be made by the user. When using
% a value of MITER = 3, the subroutine FA is
% not required, even if IMPL is not 0. For
% further information on using this option, see
% Section IV-E below.
% MITER = 4 Means the same as MITER = 1 but the A and
% Jacobian matrices are assumed to be banded.
% MITER = 5 Means the same as MITER = 2 but the A and
% Jacobian matrices are assumed to be banded.
%
% IMPL = (Input) The implicit method indicator.
% IMPL = 0 Means solving dY(I)/dT = F(Y(I),T).
% IMPL = 1 Means solving A*dY(I)/dT = F(Y(I),T), non-
% singular A (see description of FA below.)
% Only MINT = 1 or 2, and MITER = 1, 2, 3, 4,
% or 5 are allowed for this option.
% IMPL = 2,3 Means solving certain systems of hybrid
% differential/algebraic equations (see
% description of FA below.) Only MINT = 2 and
% MITER = 1, 2, 3, 4, or 5, are allowed for
% this option.
% The value of IMPL must not be changed during a problem.
%
% ML = (Input) The lower half-bandwidth in the case of a banded
% A or Jacobian matrix. (I.e., maximum(R-C) for nonzero
% A(R,C).)
%
% MU = (Input) The upper half-bandwidth in the case of a banded
% A or Jacobian matrix. (I.e., maximum(C-R).)
%
% MXORD = (Input) The maximum order desired. This is .LE. 12 for
% the Adams methods and .LE. 5 for the Gear methods. Normal
% value is 12 and 5, respectively. If MINT is 3, the
% maximum order used will be MIN(MXORD, 12) when using the
% Adams methods, and MIN(MXORD, 5) when using the Gear
% methods. MXORD must not be altered during a problem.
%
% HMAX = (Input) The maximum magnitude of the step sizemlv that will
% be used for the problem. This is useful for ensuring that
% important details are not missed. If this is not the
% case, a large value, such as the interval length, is
% suggested.
%
% WORK
% LENW = (Input)
% WORK is an array of LENW real words used
% internally for temporary storage. The user must allocate
% space for this array in the calling program by a statement
% such as
% REAL WORK(...)
% The following table gives the required minimum value for
% the length of WORK, depending on the value of IMPL and
% MITER. LENW should be set to the value used. The
% contents of WORK should not be disturbed between calls to
% SDRIV3.
%
% IMPL = 0 1 2 3
% ---------------------------------------------------------
% MITER = 0 (MXORD+4)*N Not allowed Not allowed Not allowed
% + 2*NROOT
% + 250
%
% 1,2 N*N + 2*N*N + N*N + N*(N + NDE)
% (MXORD+5)*N (MXORD+5)*N (MXORD+6)*N + (MXORD+5)*N
% + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT
% + 250 + 250 + 250 + 250
%
% 3 (MXORD+4)*N (MXORD+4)*N (MXORD+4)*N (MXORD+4)*N
% + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT
% + 250 + 250 + 250 + 250
%
% 4,5 (2*ML+MU+1) 2*(2*ML+MU+1) (2*ML+MU+1) (2*ML+MU+1)*
% *N + *N + *N + (N+NDE) +
% (MXORD+5)*N (MXORD+5)*N (MXORD+6)*N + (MXORD+5)*N
% + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT
% + 250 + 250 + 250 + 250
% ---------------------------------------------------------
%
% IWORK
% LENIW = (Input)
% IWORK is an integer array of length LENIW used internally
% for temporary storage. The user must allocate space for
% this array in the calling program by a statement such as
% INTEGER IWORK(...)
% The length of IWORK should be at least
% 50 if MITER is 0 or 3, or
% N+50 if MITER is 1, 2, 4, or 5, or MINT is 3,
% and LENIW should be set to the value used. The contents
% of IWORK should not be disturbed between calls to SDRIV3.
%
% JACOBN = A subroutine supplied by the user, if MITER is 1 or 4.
% If this is the case, the name must be declared EXTERNAL in
% the user's calling program. Given a system of N
% differential equations, it is meaningful to speak about
% the partial derivative of the I-th right hand side with
% respect to the J-th dependent variable. In general there
% are N*N such quantities. Often however the equations can
% be ordered so that the I-th differential equation only
% involves dependent variables with index near I, e.g., I+1,
% I-2. Such a system is called banded. If, for all I, the
% I-th equation depends on at most the variables
% Y(I-ML), Y(I-ML+1), ... , Y(I), Y(I+1), ... , Y(I+MU)
% then we call ML+MU+1 the bandwidth of the system. In a
% banded system many of the partial derivatives above are
% automatically zero. For the cases MITER = 1, 2, 4, and 5,
% some of these partials are needed. For the cases
% MITER = 2 and 5 the necessary derivatives are
% approximated numerically by SDRIV3, and we only ask the
% user to tell SDRIV3 the value of ML and MU if the system
% is banded. For the cases MITER = 1 and 4 the user must
% derive these partials algebraically and encode them in
% subroutine JACOBN. By computing these derivatives the
% user can often savemlv 20-30 per cent of the computing time.
% Usually, however, the accuracy is not much affected and
% most users will probably forego this option. The optional
% user-written subroutine JACOBN has the form:
% subroutine JACOBN (N, T, Y, DFDY, MATDIM, ML, MU)
% REAL Y(*), DFDY(MATDIM,*)
% .
% .
% Calculate values of DFDY
% .
% .
% end (Sample)
% Here Y is a vector of length at least N. The actual
% length of Y is determined by the user's declaration in the
% program which calls SDRIV3. Thus the dimensioning of Y in
% JACOBN, while required by FORTRAN convention, does not
% actually allocate any storage. When this subroutine is
% called, the first N components of Y are intermediate
% approximations to the solution components. The user
% should not alter these values. If the system is not
% banded (MITER=1), the partials of the I-th equation with
% respect to the J-th dependent function are to be stored in
% DFDY(I,J). Thus partials of the I-th equation are stored
% in the I-th row of DFDY. If the system is banded
% (MITER=4), then the partials of the I-th equation with
% respect to Y(J) are to be stored in DFDY(K,J), where
% K=I-J+MU+1 . Normally a return from JACOBN passes control
% back to SDRIV3. However, if the user would like to abort
% the calculation, i.e., return control to the program which
% calls SDRIV3, he should set N to zero. SDRIV3 will signal
% this by returning a value of NSTATE equal to +8(-8).
% Altering the value of N in JACOBN has no effect on the
% value of N in the call sequence of SDRIV3.
%
% FA = A subroutine supplied by the user if IMPL is not zero, and
% MITER is not 3. If so, the name must be declared EXTERNAL
% in the user's calling program. This subroutine computes
% the array A, where A*dY(I)/dT = F(Y(I),T).
% There are three cases:
%
% IMPL=1.
% subroutine FA is of the form:
% subroutine FA (N, T, Y, A, MATDIM, ML, MU, NDE)
% REAL Y(*), A(MATDIM,*)
% .
% .
% Calculate ALL values of A
% .
% .
% end (Sample)
% In this case A is assumed to be a nonsingular matrix,
% with the same structure as DFDY (see JACOBN description
% above). Programming considerations prevent complete
% generality. If MITER is 1 or 2, A is assumed to be full
% and the user must compute and store all values of
% A(I,J), I,J=1, ... ,N. If MITER is 4 or 5, A is assumed
% to be banded with lower and upper half bandwidth ML and
% MU. The left hand side of the I-th equation is a linear
% combination of dY(I-ML)/dT, dY(I-ML+1)/dT, ... ,
% dY(I)/dT, ... , dY(I+MU-1)/dT, dY(I+MU)/dT. Thus in the
% I-th equation, the coefficient of dY(J)/dT is to be
% stored in A(K,J), where K=I-J+MU+1.
% NOTE: The array A will be altered between calls to FA.
%
% IMPL=2.
% subroutine FA is of the form:
% subroutine FA (N, T, Y, A, MATDIM, ML, MU, NDE)
% REAL Y(*), A(*)
% .
% .
% Calculate non-zero values of A(1),...,A(NDE)
% .
% .
% end (Sample)
% In this case it is assumed that the system is ordered by
% the user so that the differential equations appear
% first, and the algebraic equations appear last. The
% algebraic equations must be written in the form:
% 0 = F(Y(I),T). When using this option it is up to the
% user to provide initial values for the Y(I) that satisfy
% the algebraic equations as well as possible. It is
% further assumed that A is a vector of length NDE. All
% of the components of A, which may depend on T, Y(I),
% etc., must be set by the user to non-zero values.
%
% IMPL=3.
% subroutine FA is of the form:
% subroutine FA (N, T, Y, A, MATDIM, ML, MU, NDE)
% REAL Y(*), A(MATDIM,*)
% .
% .
% Calculate ALL values of A
% .
% .
% end (Sample)
% In this case A is assumed to be a nonsingular NDE by NDE
% matrix with the same structure as DFDY (see JACOBN
% description above). Programming considerations prevent
% complete generality. If MITER is 1 or 2, A is assumed
% to be full and the user must compute and store all
% values of A(I,J), I,J=1, ... ,NDE. If MITER is 4 or 5,
% A is assumed to be banded with lower and upper half
% bandwidths ML and MU. The left hand side of the I-th
% equation is a linear combination of dY(I-ML)/dT,
% dY(I-ML+1)/dT, ... , dY(I)/dT, ... , dY(I+MU-1)/dT,
% dY(I+MU)/dT. Thus in the I-th equation, the coefficient
% of dY(J)/dT is to be stored in A(K,J), where K=I-J+MU+1.
% It is assumed that the system is ordered by the user so
% that the differential equations appear first, and the
% algebraic equations appear last. The algebraic
% equations must be written in the form 0 = F(Y(I),T).
% When using this option it is up to the user to provide
% initial values for the Y(I) that satisfy the algebraic
% equations as well as possible.
% NOTE: For IMPL = 3, the array A will be altered between
% calls to FA.
% Here Y is a vector of length at least N. The actual
% length of Y is determined by the user's declaration in the
% program which calls SDRIV3. Thus the dimensioning of Y in
% FA, while required by FORTRAN convention, does not
% actually allocate any storage. When this subroutine is
% called, the first N components of Y are intermediate
% approximations to the solution components. The user
% should not alter these values. FA is always called
% immediately after calling F, with the same values of T
% and Y. Normally a return from FA passes control back to
% SDRIV3. However, if the user would like to abort the
% calculation, i.e., return control to the program which
% calls SDRIV3, he should set N to zero. SDRIV3 will signal
% this by returning a value of NSTATE equal to +9(-9).
% Altering the value of N in FA has no effect on the value
% of N in the call sequence of SDRIV3.
%
% NDE = (Input) The number of differential equations. This is
% required only for IMPL = 2 or 3, with NDE .LT. N.
%
% MXSTEP = (Input) The maximum number of internal steps allowed on
% one call to SDRIV3.
%
% G = A real FORTRAN function supplied by the user
% if NROOT is not 0. In this case, the name must be
% declared EXTERNAL in the user's calling program. G is
% repeatedly called with different values of IROOT to obtain
% the value of each of the NROOT equations for which a root
% is desired. G is of the form:
% REAL function G (N, T, Y, IROOT)
% REAL Y(*)
% GO TO (10, ...), IROOT
% 10 G = ...
% .
% .
% end (Sample)
% Here, Y is a vector of length at least N, whose first N
% components are the solution components at the point T.
% The user should not alter these values. The actual length
% of Y is determined by the user's declaration in the
% program which calls SDRIV3. Thus the dimensioning of Y in
% G, while required by FORTRAN convention, does not actually
% allocate any storage. Normally a return from G passes
% control back to SDRIV3. However, if the user would like
% to abort the calculation, i.e., return control to the
% program which calls SDRIV3, he should set N to zero.
% SDRIV3 will signal this by returning a value of NSTATE
% equal to +7(-7). In this case, the index of the equation
% being evaluated is stored in the sixth element of IWORK.
% Altering the value of N in G has no effect on the value of
% N in the call sequence of SDRIV3.
%
% USERS = A subroutine supplied by the user, if MITER is 3.
% If this is the case, the name must be declared EXTERNAL in
% the user's calling program. The routine USERS is called
% by SDRIV3 when certain linear systems must be solved. The
% user may choose any method to form, store and solve these
% systems in order to obtain the solution result that is
% returned to SDRIV3. In particular, this allows sparse
% matrix methods to be used. The call sequence for this
% routine is:
%
% subroutine USERS (Y, YH, YWT, SAVE1, SAVE2, T, H, EL,
% 8 IMPL, N, NDE, IFLAG)
% REAL Y(*), YH(*), YWT(*), SAVE1(*),
% 8 SAVE2(*), T, H, EL
%
% The input variable IFLAG indicates what action is to be
% taken. Subroutine USERS should perform the following
% operations, depending on the value of IFLAG and IMPL.
%
% IFLAG = 0
% IMPL = 0. USERS is not called.
% IMPL = 1, 2 or 3. Solve the system A*X = SAVE2,
% returning the result in SAVE2. The array SAVE1 can
% be used as a work array. For IMPL = 1, there are N
% components to the system, and for IMPL = 2 or 3,
% there are NDE components to the system.
%
% IFLAG = 1
% IMPL = 0. Compute, decompose and store the matrix
% (I - H*EL*J), where I is the identity matrix and J
% is the Jacobian matrix of the right hand side. The
% array SAVE1 can be used as a work array.
% IMPL = 1, 2 or 3. Compute, decompose and store the
% matrix (A - H*EL*J). The array SAVE1 can be used as
% a work array.
%
% IFLAG = 2
% IMPL = 0. Solve the system
% (I - H*EL*J)*X = H*SAVE2 - YH - SAVE1,
% returning the result in SAVE2.
% IMPL = 1, 2 or 3. Solve the system
% (A - H*EL*J)*X = H*SAVE2 - A*(YH + SAVE1)
% returning the result in SAVE2.
% The array SAVE1 should not be altered.
% If IFLAG is 0 and IMPL is 1 or 2 and the matrix A is
% singular, or if IFLAG is 1 and one of the matrices
% (I - H*EL*J), (A - H*EL*J) is singular, the INTEGER
% variable IFLAG is to be set to -1 before RETURNing.
% Normally a return from USERS passes control back to
% SDRIV3. However, if the user would like to abort the
% calculation, i.e., return control to the program which
% calls SDRIV3, he should set N to zero. SDRIV3 will signal
% this by returning a value of NSTATE equal to +10(-10).
% Altering the value of N in USERS has no effect on the
% value of N in the call sequence of SDRIV3.
%
% IERFLG = An error flag. The error number associated with a
% diagnostic message (see Section III-A below) is the same
% as the corresponding value of IERFLG. The meaning of
% IERFLG:
% 0 The routine completed successfully. (No message is
% issued.)
% 3 (Warning) The number of steps required to reach TOUT
% exceeds MXSTEP.
% 4 (Warning) The value of EPS is too small.
% 11 (Warning) For NTASK = 2 or 3, T is beyond TOUT.
% The solution was obtained by interpolation.
% 15 (Warning) The integration step sizemlv is below the
% roundoff level of T. (The program issues this
% message as a warning but does not return control to
% the user.)
% 22 (Recoverable) N is not positive.
% 23 (Recoverable) MINT is less than 1 or greater than 3 .
% 24 (Recoverable) MITER is less than 0 or greater than
% 5 .
% 25 (Recoverable) IMPL is less than 0 or greater than 3 .
% 26 (Recoverable) The value of NSTATE is less than 1 or
% greater than 12 .
% 27 (Recoverable) EPS is less than zero.
% 28 (Recoverable) MXORD is not positive.
% 29 (Recoverable) For MINT = 3, either MITER = 0 or 3, or
% IMPL = 0 .
% 30 (Recoverable) For MITER = 0, IMPL is not 0 .
% 31 (Recoverable) For MINT = 1, IMPL is 2 or 3 .
% 32 (Recoverable) Insufficient storage has been allocated
% for the WORK array.
% 33 (Recoverable) Insufficient storage has been allocated
% for the IWORK array.
% 41 (Recoverable) The integration step sizemlv has gone
% to zero.
% 42 (Recoverable) The integration step sizemlv has been
% reduced about 50 times without advancing the
% solution. The problem setup may not be correct.
% 43 (Recoverable) For IMPL greater than 0, the matrix A
% is singular.
% 999 (Fatal) The value of NSTATE is 12 .
%
% III. OTHER COMMUNICATION TO THE USER ..............................
%
% A. The solver communicates to the user through the parameters
% above. In addition it writes diagnostic messages through the
% standard error handling program XERMSG. A complete description
% of XERMSG is given in 'Guide to the SLATEC Common Mathematical
% Library' by Kirby W. Fong et al.. At installations which do not
% have this error handling package the short but serviceable
% routine, XERMSG, available with this package, can be used. That
% program uses the file named OUTPUT to transmit messages.
%
% B. The first three elements of WORK and the first five elements of
% IWORK will contain the following statistical data:
% AVGH The average step sizemlv used.
% HUSED The step sizemlv last used (successfully).
% AVGORD The average order used.
% IMXERR The index of the element of the solution vector that
% contributed most to the last error test.
% NQUSED The order last used (successfully).
% NSTEP The number of steps taken since last initialization.
% NFE The number of evaluations of the right hand side.
% NJE The number of evaluations of the Jacobian matrix.
%
% IV. REMARKS .......................................................
%
% A. Other routines used:
% SDNTP, SDZRO, SDSTP, SDNTL, SDPST, SDCOR, SDCST,
% SDPSC, and SDSCL;
% SGEFA, SGESL, SGBFA, SGBSL, and SNRM2 (from LINPACK)
% R1MACH (from the Bell Laboratories Machine Constants Package)
% XERMSG (from the SLATEC Common Math Library)
% The last seven routines above, not having been written by the
% present authors, are not explicitly part of this package.
%
% B. On any return from SDRIV3 all information necessary to continue
% the calculation is contained in the call sequence parameters,
% including the work arrays. Thus it is possible to suspend one
% problem, integrate another, and then return to the first.
%
% C. If this package is to be used in an overlay situation, the user
% must declare in the primary overlay the variables in the call
% sequence to SDRIV3.
%
% D. Changing parameters during an integration.
% The value of NROOT, EPS, EWT, IERROR, MINT, MITER, or HMAX may
% be altered by the user between calls to SDRIV3. For example, if
% too much accuracy has been requested (the program returns with
% NSTATE = 4 and an increased value of EPS) the user may wish to
% increase EPS further. In general, prudence is necessary when
% making changes in parameters since such changes are not
% implemented until the next integration step, which is not
% necessarily the next call to SDRIV3. This can happen if the
% program has already integrated to a point which is beyond the
% new point TOUT.
%
% E. As the price for complete control of matrix algebra, the SDRIV3
% USERS option puts all responsibility for Jacobian matrix
% evaluation on the user. It is often useful to approximate
% numerically all or part of the Jacobian matrix. However this
% must be done carefully. The FORTRAN sequence below illustrates
% the method we recommend. It can be inserted directly into
% subroutine USERS to approximate Jacobian elements in rows I1
% to I2 and columns J1 to J2.
% REAL DFDY(N,N), EPSJ, H, R, R1MACH,
% 8 SAVE1(N), SAVE2(N), T, UROUND, Y(N), YJ, YWT(N)
% UROUND = R1MACH(4)
% EPSJ = SQRT(UROUND)
% DO 30 J = J1,J2
% R = EPSJ*MAX(ABS(YWT(J)), ABS(Y(J)))
% IF (R .EQ. 0.0E0) R = YWT(J)
% YJ = Y(J)
% Y(J) = Y(J) + R
% CALL F (N, T, Y, SAVE1)
% IF (N .EQ. 0) RETURN
% Y(J) = YJ
% DO 20 I = I1,I2
% 20 DFDY(I,J) = (SAVE1(I) - SAVE2(I))/R
% 30 CONTINUE
% Many problems give rise to structured sparse Jacobians, e.g.,
% block banded. It is possible to approximate them with fewer
% function evaluations than the above procedure uses; see Curtis,
% Powell and Reid, J. Inst. Maths Applics, (1974), Vol. 13,
% pp. 117-119.
%
% F. When any of the routines JACOBN, FA, G, or USERS, is not
% required, difficulties associated with unsatisfied externals can
% be avoided by using the name of the routine which calculates the
% right hand side of the differential equations in place of the
% corresponding name in the call sequence of SDRIV3.
%
%***REFERENCES C. W. Gear, Numerical Initial Value Problems in
% Ordinary Differential Equations, Prentice-Hall, 1971.
%***ROUTINES CALLED R1MACH, SDNTP, SDSTP, SDZRO, SGBFA, SGBSL, SGEFA,
% SGESL, SNRM2, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790601 DATE WRITTEN
% 900329 Initial submission to SLATEC.
%***end PROLOGUE SDRIV3
persistent ae big convrg glast gnow gt h hsignmlv hused i ia iavgh iavgrd icnvrg idfdy iel ifac iflag ignow ih ihmax ihold ihsign ihused ii ijroot ijstpl ijtask imach1 imach4 imnt imntld imtr imtrld imtrsv imxerr imxord imxrds indmxr indprt indpvt indtrt infe info inje inq inquse inroot inrtld instep intgr1 intgr2 inwait irc irmax iroot isave1 isave2 it itout itq itrend itroot iyh iywt j jstate jtroot lenchk liwchk matdim maxord ndecom npar nround nstepl re rl1 rl2 sizemlv summlv tlast troot uround ;
if isempty(ae), ae=0; end;
if isempty(big), big=0; end;
ewt_shape=size(ewt);ewt=reshape(ewt,1,[]);
if isempty(glast), glast=0; end;
if isempty(gnow), gnow=0; end;
if isempty(h), h=0; end;
if isempty(hsignmlv), hsignmlv=0; end;
if isempty(hused), hused=0; end;
if isempty(re), re=0; end;
if isempty(sizemlv), sizemlv=0; end;
if isempty(summlv), summlv=0; end;
if isempty(tlast), tlast=0; end;
if isempty(troot), troot=0; end;
if isempty(uround), uround=0; end;
work_shape=size(work);work=reshape(work,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
if isempty(i), i=0; end;
if isempty(ia), ia=0; end;
if isempty(idfdy), idfdy=0; end;
if isempty(ifac), ifac=0; end;
if isempty(iflag), iflag=0; end;
if isempty(ignow), ignow=0; end;
if isempty(imxerr), imxerr=0; end;
if isempty(info), info=0; end;
if isempty(iroot), iroot=0; end;
if isempty(isave1), isave1=0; end;
if isempty(isave2), isave2=0; end;
if isempty(itroot), itroot=0; end;
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
if isempty(iywt), iywt=0; end;
if isempty(j), j=0; end;
if isempty(jstate), jstate=0; end;
if isempty(jtroot), jtroot=0; end;
if isempty(lenchk), lenchk=0; end;
if isempty(liwchk), liwchk=0; end;
if isempty(matdim), matdim=0; end;
if isempty(maxord), maxord=0; end;
if isempty(ndecom), ndecom=0; end;
if isempty(npar), npar=0; end;
if isempty(nstepl), nstepl=0; end;
if isempty(ii), ii=0; end;
if isempty(gt), gt=zeros(1,7); end;
if isempty(convrg), convrg=false; end;
if isempty(intgr1), intgr1=repmat(' ',1,8); end;
if isempty(intgr2), intgr2=repmat(' ',1,8); end;
if isempty(rl1), rl1=repmat(' ',1,16); end;
if isempty(rl2), rl2=repmat(' ',1,16); end;
if isempty(nround), nround=20.0e0 ; end;
if isempty(iavgh), iavgh=1; end;
if isempty(ihused), ihused=2; end;
if isempty(iavgrd), iavgrd=3; end;
if isempty(iel), iel=4; end;
if isempty(ih), ih=160; end;
if isempty(ihmax), ihmax=161; end;
if isempty(ihold), ihold=162; end;
if isempty(ihsign), ihsign=163; end;
if isempty(irc), irc=164; end;
if isempty(irmax), irmax=165; end;
if isempty(it), it=166; end;
if isempty(itout), itout=167; end;
if isempty(itq), itq=168; end;
if isempty(itrend), itrend=204; end;
if isempty(imach1), imach1=205; end;
if isempty(imach4), imach4=206; end;
if isempty(iyh), iyh=251; end;
if isempty(indmxr), indmxr=1; end;
if isempty(inquse), inquse=2; end;
if isempty(instep), instep=3; end;
if isempty(infe), infe=4; end;
if isempty(inje), inje=5; end;
if isempty(inroot), inroot=6; end;
if isempty(icnvrg), icnvrg=7; end;
if isempty(ijroot), ijroot=8; end;
if isempty(ijtask), ijtask=9; end;
if isempty(imntld), imntld=10; end;
if isempty(imtrld), imtrld=11; end;
if isempty(inq), inq=12; end;
if isempty(inrtld), inrtld=13; end;
if isempty(indtrt), indtrt=14; end;
if isempty(inwait), inwait=15; end;
if isempty(imnt), imnt=16; end;
if isempty(imtrsv), imtrsv=17; end;
if isempty(imtr), imtr=18; end;
if isempty(imxrds), imxrds=19; end;
if isempty(imxord), imxord=20; end;
if isempty(indprt), indprt=21; end;
if isempty(ijstpl), ijstpl=22; end;
if isempty(indpvt), indpvt=51 ; end;
%***FIRST EXECUTABLE STATEMENT SDRIV3
if( nstate==12 )
ierflg = 999;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3','Illegal input. The value of NSTATE is 12 .',ierflg,2);
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( nstate<1 || nstate>12 ) ;
intgr1=sprintf(['%8i'], nstate);
ierflg = 26;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Improper value for NSTATE(= ',[intgr1,').']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
npar = fix(n);
if( eps<0.0e0 )
rl1=sprintf(['%16.8f'], eps);
ierflg = 27;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. EPS, ',[rl1,', is negative.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( n<=0 )
intgr1=sprintf(['%8i'], n);
ierflg = 22;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Number of equations, ',[intgr1,', is not positive.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( mxord<=0 )
intgr1=sprintf(['%8i'], mxord);
ierflg = 28;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Maximum order, ',[intgr1,', is not positive.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( mint<1 || mint>3 )
intgr1=sprintf(['%8i'], mint);
ierflg = 23;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Improper value for the integration method ',['flag, ',[intgr1,' .']]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( miter<0 || miter>5 ) ;
intgr1=sprintf(['%8i'], miter);
ierflg = 24;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Improper value for MITER(= ',[intgr1,').']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( impl<0 || impl>3 ) ;
intgr1=sprintf(['%8i'], impl);
ierflg = 25;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Improper value for IMPL(= ',[intgr1,').']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( mint==3 &&(miter==0 || miter==3 || impl~=0) );
intgr1=sprintf(['%8i'], miter);
intgr2=sprintf(['%8i'], impl);
ierflg = 29;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. For MINT = 3, the value of MITER, ',[intgr1,[', and/or IMPL, ',[intgr2,', is not allowed.']]]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif((impl>=1 && impl<=3) && miter==0 ) ;
intgr1=sprintf(['%8i'], impl);
ierflg = 30;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. For MITER = 0, the value of IMPL, ',[intgr1,', is not allowed.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif((impl==2 || impl==3) && mint==1 ) ;
intgr1=sprintf(['%8i'], impl);
ierflg = 31;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. For MINT = 1, the value of IMPL, ',[intgr1,', is not allowed.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( miter==0 || miter==3 )
liwchk = fix(indpvt - 1);
elseif( miter==1 || miter==2 || miter==4 || miter==5 ) ;
liwchk = fix(indpvt + n - 1);
end;
if( leniw<liwchk )
intgr1=sprintf(['%8i'], liwchk);
ierflg = 33;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Insufficient storage allocated for the ',['IWORK array. Based on the value of the input parameters ',['involved, the required storage is ',[intgr1,' .']]]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
% Allocate the WORK array
% IYH is the index of YH in WORK
if( mint==1 || mint==3 )
maxord = fix(min(mxord,12));
elseif( mint==2 ) ;
maxord = fix(min(mxord,5));
end;
idfdy = fix(iyh +(maxord+1).*n);
% IDFDY is the index of DFDY
%
if( miter==0 || miter==3 )
iywt = fix(idfdy);
elseif( miter==1 || miter==2 ) ;
iywt = fix(idfdy + n.*n);
elseif( miter==4 || miter==5 ) ;
iywt = fix(idfdy +(2.*ml+mu+1).*n);
end;
% IYWT is the index of YWT
isave1 = fix(iywt + n);
% ISAVE1 is the index of SAVE1
isave2 = fix(isave1 + n);
% ISAVE2 is the index of SAVE2
ignow = fix(isave2 + n);
% IGNOW is the index of GNOW
itroot = fix(ignow + nroot);
% ITROOT is the index of TROOT
ifac = fix(itroot + nroot);
% IFAC is the index of FAC
if( miter==2 || miter==5 || mint==3 )
ia = fix(ifac + n);
else;
ia = fix(ifac);
end;
% IA is the index of A
if( impl==0 || miter==3 )
lenchk = fix(ia - 1);
elseif( impl==1 &&(miter==1 || miter==2) ) ;
lenchk = fix(ia - 1 + n.*n);
elseif( impl==1 &&(miter==4 || miter==5) ) ;
lenchk = fix(ia - 1 +(2.*ml+mu+1).*n);
elseif( impl==2 && miter~=3 ) ;
lenchk = fix(ia - 1 + n);
elseif( impl==3 &&(miter==1 || miter==2) ) ;
lenchk = fix(ia - 1 + n.*nde);
elseif( impl==3 &&(miter==4 || miter==5) ) ;
lenchk = fix(ia - 1 +(2.*ml+mu+1).*nde);
end;
if( lenw<lenchk )
intgr1=sprintf(['%8i'], lenchk);
ierflg = 32;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['Illegal input. Insufficient storage allocated for the ',['WORK array. Based on the value of the input parameters ',['involved, the required storage is ',[intgr1,' .']]]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
if( miter==0 || miter==3 )
matdim = 1;
elseif( miter==1 || miter==2 ) ;
matdim = fix(n);
elseif( miter==4 || miter==5 ) ;
matdim = fix(2.*ml + mu + 1);
end;
if( impl==0 || impl==1 )
ndecom = fix(n);
elseif( impl==2 || impl==3 ) ;
ndecom = fix(nde);
end;
gt(:)=0;
while (1);
if( nstate==1 )
% Initialize parameters
if( mint==1 || mint==3 )
iwork(imxord) = fix(min(mxord,12));
elseif( mint==2 ) ;
iwork(imxord) = fix(min(mxord,5));
end;
iwork(imxrds) = fix(mxord);
if( mint==1 || mint==2 )
iwork(imnt) = fix(mint);
iwork(imtr) = fix(miter);
iwork(imntld) = fix(mint);
iwork(imtrld) = fix(miter);
elseif( mint==3 ) ;
iwork(imnt) = 1;
iwork(imtr) = 0;
iwork(imntld) = fix(iwork(imnt));
iwork(imtrld) = fix(iwork(imtr));
iwork(imtrsv) = fix(miter);
end;
work(ihmax) = hmax;
[uround ]=r1mach(4);
work(imach4) = uround;
[work(imach1) ]=r1mach(1);
if( nroot~=0 )
re = uround;
ae = work(imach1);
end;
h =(tout-t).*(1.0e0-4.0e0.*uround);
h = (abs(min(abs(h),hmax)).*sign(h));
work(ih) = h;
hsignmlv = (abs(1.0e0).*sign(h));
work(ihsign) = hsignmlv;
iwork(ijtask) = 0;
work(iavgh) = 0.0e0;
work(ihused) = 0.0e0;
work(iavgrd) = 0.0e0;
iwork(indmxr) = 0;
iwork(inquse) = 0;
iwork(instep) = 0;
iwork(ijstpl) = 0;
iwork(infe) = 0;
iwork(inje) = 0;
iwork(inroot) = 0;
work(it) = t;
iwork(icnvrg) = 0;
iwork(indprt) = 0;
% Set initial conditions
for i = 1 : n;
work(i+iyh-1) = y(i);
end; i = fix(n+1);
if( t==tout )
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
break;
else;
uround = work(imach4);
if( nroot~=0 )
re = uround;
ae = work(imach1);
end;
end;
% On a continuation, check
% that output points have
% been or will be overtaken.
if( iwork(icnvrg)==1 )
convrg = true;
else;
convrg = false;
end;
t = work(it);
h = work(ih);
hsignmlv = work(ihsign);
if( iwork(ijtask)~=0 )
%
% IWORK(IJROOT) flags unreported
% roots, and is set to the value of
% NTASK when a root was last selected.
% It is set to zero when all roots
% have been reported. IWORK(INROOT)
% contains the index and WORK(ITOUT)
% contains the value of the root last
% selected to be reported.
% IWORK(INRTLD) contains the value of
% NROOT and IWORK(INDTRT) contains
% the value of ITROOT when the array
% of roots was last calculated.
if( nroot~=0 )
if( iwork(ijroot)>0 )
% TOUT has just been reported.
% If TROOT .LE. TOUT, report TROOT.
if( nstate==5 )
troot = t;
iroot = 0;
for i = 1 : iwork(inrtld);
jtroot = fix(i + iwork(indtrt) - 1);
if( work(jtroot).*hsignmlv<=troot.*hsignmlv )
%
% Check for multiple roots.
%
if( work(jtroot)==work(itout) && i>iwork(inroot) )
iroot = fix(i);
troot = work(jtroot);
break;
end;
if( work(jtroot).*hsignmlv>work(itout).*hsignmlv )
iroot = fix(i);
troot = work(jtroot);
end;
end;
end;
iwork(inroot) = fix(iroot);
work(itout) = troot;
iwork(ijroot) = fix(ntask);
if( ntask==1 )
if( iroot==0 )
iwork(ijroot) = 0;
elseif( tout.*hsignmlv>=troot.*hsignmlv ) ;
[h,dumvar2,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y);
nstate = 5;
t = troot;
ierflg = 0;
gt(5)=1;
break;
end;
elseif( ntask==2 || ntask==3 ) ;
%
% If there are no more roots, or the
% user has altered TOUT to be less
% than a root, set IJROOT to zero.
%
if( iroot==0 ||(tout.*hsignmlv<troot.*hsignmlv) )
iwork(ijroot) = 0;
else;
[h,dumvar2,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y);
nstate = 5;
ierflg = 0;
t = troot;
gt(5)=1;
break;
end;
end;
elseif( tout.*hsignmlv>=work(itout).*hsignmlv ) ;
troot = work(itout);
[h,dumvar2,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y);
t = troot;
nstate = 5;
ierflg = 0;
gt(5)=1;
break;
% A root has just been reported.
% Select the next root.
end;
end;
end;
%
if( ntask==1 )
nstate = 2;
if( t.*hsignmlv>=tout.*hsignmlv )
[h,dumvar2,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y);
t = tout;
ierflg = 0;
gt(5)=1;
break;
end;
elseif( ntask==2 ) ;
% Check if TOUT has
% been reset .LT. T
if( t.*hsignmlv>tout.*hsignmlv )
rl1=sprintf(['%16.8f'], t);
rl2=sprintf(['%16.8f'], tout);
ierflg = 11;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['While integrating exactly to TOUT, T, ',[rl1,[', was beyond TOUT, ',[rl2,[' . Solution obtained by ','interpolation.']]]]],ierflg,0);
nstate = 11;
[h,dumvar2,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y);
t = tout;
gt(5)=1;
break;
end;
% Determine if TOUT has been overtaken
%
if( abs(tout-t)<=nround.*uround.*max(abs(t),abs(tout)) )
t = tout;
nstate = 2;
ierflg = 0;
gt(4)=1;
break;
end;
% If there are no more roots
% to report, report T.
if( nstate==5 )
nstate = 2;
ierflg = 0;
gt(4)=1;
break;
end;
nstate = 2;
% See if TOUT will
% be overtaken.
if((t+h).*hsignmlv>tout.*hsignmlv )
h = tout - t;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = h.*(1.0e0-4.0e0.*uround);
end;
work(ih) = h;
if( h==0.0e0 )
gt(6)=1;
break;
end;
iwork(ijtask) = -1;
end;
elseif( ntask==3 ) ;
nstate = 2;
if( t.*hsignmlv>tout.*hsignmlv )
rl1=sprintf(['%16.8f'], t);
rl2=sprintf(['%16.8f'], tout);
ierflg = 11;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['While integrating exactly to TOUT, T, ',[rl1,[', was beyond TOUT, ',[rl2,[' . Solution obtained by ','interpolation.']]]]],ierflg,0);
nstate = 11;
[h,dumvar2,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y);
t = tout;
gt(5)=1;
break;
end;
if( abs(tout-t)<=nround.*uround.*max(abs(t),abs(tout)) )
t = tout;
ierflg = 0;
gt(4)=1;
break;
end;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = tout - t;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = h.*(1.0e0-4.0e0.*uround);
end;
work(ih) = h;
gt(6)=1;
break;
iwork(ijtask) = -1;
end;
end;
% Implement changes in MINT, MITER, and/or HMAX.
%
if((mint~=iwork(imntld) || miter~=iwork(imtrld)) &&mint~=3 && iwork(imntld)~=3 )
iwork(ijtask) = -1;
end;
if( hmax~=work(ihmax) )
h = (abs(min(abs(h),hmax)).*sign(h));
if( h~=work(ih) )
iwork(ijtask) = -1;
work(ih) = h;
end;
work(ihmax) = hmax;
end;
end;
%
break;
end;
while (1);
if(gt(6)==0)
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
nstepl = fix(iwork(instep));
for i = 1 : n;
y(i) = work(i+iyh-1);
end; i = fix(n+1);
if( nroot~=0 )
for i = 1 : nroot;
work(i+ignow-1) = g(npar,t,y,i);
if( npar==0 )
iwork(inroot) = fix(i);
nstate = 7;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
end; i = fix(nroot+1);
end;
if( ierror==1 )
for i = 1 : n;
work(i+iywt-1) = 1.0e0;
end; i = fix(n+1);
gt(3)=1;
continue;
elseif( ierror==5 ) ;
for i = 1 : n;
work(i+iywt-1) = ewt(i);
end; i = fix(n+1);
gt(3)=1;
continue;
end;
% Reset YWT array. Looping point.
end;
gt(2)=0;
if( ierror==2 )
for i = 1 : n;
if( y(i)==0.0e0 )
if( iwork(ijtask)==0 )
[npar,t,y,work(isave2:end)]=f(npar,t,y,work(isave2:end));
if( npar==0 )
nstate = 6;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
iwork(infe) = fix(iwork(infe) + 1);
if( miter==3 && impl~=0 )
iflag = 0;
[y,dumvar2,dumvar3,dumvar4,dumvar5,t,h,dumvar8,impl,npar,ndecom,iflag]=users(y,work(iyh:end),work(iywt:end),work(isave1:end),work(isave2:end),t,h,work(iel:end),impl,npar,ndecom,iflag); dumvar2i=find((work(iyh:end))~=(dumvar2));dumvar3i=find((work(iywt:end))~=(dumvar3));dumvar4i=find((work(isave1:end))~=(dumvar4));dumvar5i=find((work(isave2:end))~=(dumvar5));dumvar8i=find((work(iel:end))~=(dumvar8)); work(iyh-1+dumvar2i)=dumvar2(dumvar2i); work(iywt-1+dumvar3i)=dumvar3(dumvar3i); work(isave1-1+dumvar4i)=dumvar4(dumvar4i); work(isave2-1+dumvar5i)=dumvar5(dumvar5i); work(iel-1+dumvar8i)=dumvar8(dumvar8i);
if( iflag==-1 )
gt(7)=1;
break;
end;
if( npar==0 )
nstate = 10;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
elseif( impl==1 ) ;
if( miter==1 || miter==2 )
[npar,t,y,work(ia:end),matdim,ml,mu,ndecom]=fa(npar,t,y,work(ia:end),matdim,ml,mu,ndecom);
if( npar==0 )
nstate = 9;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
[work(sub2ind(size(work),max(ia,1)):end),matdim,n,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info]=sgefa(work(sub2ind(size(work),max(ia,1)):end),matdim,n,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info);
if( info~=0 )
gt(7)=1;
break;
end;
[dumvar1,matdim,n,iwork(sub2ind(size(iwork),max(indpvt,1)):end),dumvar5]=sgesl(work(sub2ind(size(work),max(ia,1)):end),matdim,n,iwork(sub2ind(size(iwork),max(indpvt,1)):end),work(sub2ind(size(work),max(isave2,1)):end),0); dumvar1i=find((work(sub2ind(size(work),max(ia,1)):end))~=(dumvar1));dumvar5i=find((work(sub2ind(size(work),max(isave2,1)):end))~=(dumvar5)); work(ia-1+dumvar1i)=dumvar1(dumvar1i); work(isave2-1+dumvar5i)=dumvar5(dumvar5i);
elseif( miter==4 || miter==5 ) ;
[npar,t,y,work(ia+ml:end),matdim,ml,mu,ndecom]=fa(npar,t,y,work(ia+ml:end),matdim,ml,mu,ndecom);
if( npar==0 )
nstate = 9;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
[work(sub2ind(size(work),max(ia,1)):end),matdim,n,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info]=sgbfa(work(sub2ind(size(work),max(ia,1)):end),matdim,n,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info);
if( info~=0 )
gt(7)=1;
break;
end;
[dumvar1,matdim,n,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),dumvar7]=sgbsl(work(sub2ind(size(work),max(ia,1)):end),matdim,n,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),work(sub2ind(size(work),max(isave2,1)):end),0); dumvar1i=find((work(sub2ind(size(work),max(ia,1)):end))~=(dumvar1));dumvar7i=find((work(sub2ind(size(work),max(isave2,1)):end))~=(dumvar7)); work(ia-1+dumvar1i)=dumvar1(dumvar1i); work(isave2-1+dumvar7i)=dumvar7(dumvar7i);
end;
elseif( impl==2 ) ;
[npar,t,y,work(ia:end),matdim,ml,mu,ndecom]=fa(npar,t,y,work(ia:end),matdim,ml,mu,ndecom);
if( npar==0 )
nstate = 9;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
for ii = 1 : ndecom;
if( work(ii+ia-1)==0.0e0 )
gt(7)=1;
break;
end;
work(ii+isave2-1) = work(ii+isave2-1)./work(ii+ia-1);
end;
if(gt(7)==1)
break;
end;
elseif( impl==3 ) ;
if( miter==1 || miter==2 )
[npar,t,y,work(ia:end),matdim,ml,mu,ndecom]=fa(npar,t,y,work(ia:end),matdim,ml,mu,ndecom);
if( npar==0 )
nstate = 9;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
[work(sub2ind(size(work),max(ia,1)):end),matdim,nde,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info]=sgefa(work(sub2ind(size(work),max(ia,1)):end),matdim,nde,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info);
if( info~=0 )
gt(7)=1;
break;
end;
[dumvar1,matdim,nde,iwork(sub2ind(size(iwork),max(indpvt,1)):end),dumvar5]=sgesl(work(sub2ind(size(work),max(ia,1)):end),matdim,nde,iwork(sub2ind(size(iwork),max(indpvt,1)):end),work(sub2ind(size(work),max(isave2,1)):end),0); dumvar1i=find((work(sub2ind(size(work),max(ia,1)):end))~=(dumvar1));dumvar5i=find((work(sub2ind(size(work),max(isave2,1)):end))~=(dumvar5)); work(ia-1+dumvar1i)=dumvar1(dumvar1i); work(isave2-1+dumvar5i)=dumvar5(dumvar5i);
elseif( miter==4 || miter==5 ) ;
[npar,t,y,work(ia+ml:end),matdim,ml,mu,ndecom]=fa(npar,t,y,work(ia+ml:end),matdim,ml,mu,ndecom);
if( npar==0 )
nstate = 9;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
[work(sub2ind(size(work),max(ia,1)):end),matdim,nde,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info]=sgbfa(work(sub2ind(size(work),max(ia,1)):end),matdim,nde,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),info);
if( info~=0 )
gt(7)=1;
break;
end;
[dumvar1,matdim,nde,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),dumvar7]=sgbsl(work(sub2ind(size(work),max(ia,1)):end),matdim,nde,ml,mu,iwork(sub2ind(size(iwork),max(indpvt,1)):end),work(sub2ind(size(work),max(isave2,1)):end),0); dumvar1i=find((work(sub2ind(size(work),max(ia,1)):end))~=(dumvar1));dumvar7i=find((work(sub2ind(size(work),max(isave2,1)):end))~=(dumvar7)); work(ia-1+dumvar1i)=dumvar1(dumvar1i); work(isave2-1+dumvar7i)=dumvar7(dumvar7i);
end;
end;
end;
for j = i : n;
if( y(j)~=0.0e0 )
work(j+iywt-1) = abs(y(j));
elseif( iwork(ijtask)==0 ) ;
work(j+iywt-1) = abs(h.*work(j+isave2-1));
else;
work(j+iywt-1) = abs(work(j+iyh+n-1));
end;
if( work(j+iywt-1)==0.0e0 )
work(j+iywt-1) = uround;
end;
end; j = fix(n+1);
break;
else;
work(i+iywt-1) = abs(y(i));
end;
end;
if(gt(7)==1)
break;
end;
elseif( ierror==3 ) ;
for i = 1 : n;
work(i+iywt-1) = max(ewt(1),abs(y(i)));
end; i = fix(n+1);
elseif( ierror==4 ) ;
for i = 1 : n;
work(i+iywt-1) = max(ewt(i),abs(y(i)));
end; i = fix(n+1);
end;
%
end;
gt(3)=0;
for i = 1 : n;
work(i+isave2-1) = y(i)./work(i+iywt-1);
end; i = fix(n+1);
summlv = snrm2(n,work(sub2ind(size(work),max(isave2,1)):end),1)./sqrt(real(n));
summlv = max(1.0e0,summlv);
if( eps<summlv.*uround )
eps = summlv.*uround.*(1.0e0+10.0e0.*uround);
rl1=sprintf(['%16.8f'], t);
rl2=sprintf(['%16.8f'], eps);
ierflg = 4;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,[', the requested accuracy, EPS, was not ',['obtainable with the machine precision. EPS has been ',['increased to ',[rl2,' .']]]]]],ierflg,0);
nstate = 4;
gt(4)=1;
continue;
end;
if( abs(h)>=uround.*abs(t) )
iwork(indprt) = 0;
elseif( iwork(indprt)==0 ) ;
rl1=sprintf(['%16.8f'], t);
rl2=sprintf(['%16.8f'], h);
ierflg = 15;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,[', the step sizemlv, ',[rl2,[', is smaller ',['than the roundoff level of T. This may occur if there is ',['an abrupt change in the right hand side of the ','differential equations.']]]]]]],ierflg,0);
iwork(indprt) = 1;
end;
if( ntask~=2 )
if((iwork(instep)-nstepl)==mxstep )
rl1=sprintf(['%16.8f'], t);
intgr1=sprintf(['%8i'], mxstep);
rl2=sprintf(['%16.8f'], tout);
ierflg = 3;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,[', ',[intgr1,[' steps have been taken ',['without reaching TOUT, ',[rl2,' .']]]]]]],ierflg,0);
nstate = 3;
gt(4)=1;
continue;
end;
end;
%
% CALL SDSTP (EPS, F, FA, HMAX, IMPL, IERROR, JACOBN, MATDIM,
% 8 MAXORD, MINT, MITER, ML, MU, N, NDE, YWT, UROUND,
% 8 USERS, AVGH, AVGORD, H, HUSED, JTASK, MNTOLD, MTROLD,
% 8 NFE, NJE, NQUSED, NSTEP, T, Y, YH, A, CONVRG,
% 8 DFDY, EL, FAC, HOLD, IPVT, JSTATE, JSTEPL, NQ, NWAIT,
% 8 RC, RMAX, SAVE1, SAVE2, TQ, TREND, ISWFLG, MTRSV,
% 8 MXRDSV)
%
[eps,f,fa,dumvar4,impl,ierror,jacobn,matdim,dumvar9,dumvar10,dumvar11,ml,mu,npar,ndecom,dumvar16,uround,users,dumvar19,dumvar20,dumvar21,hused,dumvar23,dumvar24,dumvar25,dumvar26,dumvar27,dumvar28,dumvar29,dumvar30,y,dumvar32,dumvar33,convrg,dumvar35,dumvar36,dumvar37,dumvar38,dumvar39,jstate,dumvar41,dumvar42,dumvar43,dumvar44,dumvar45,dumvar46,dumvar47,dumvar48,dumvar49,mint,dumvar51,dumvar52]=sdstp(eps,f,fa,work(ihmax),impl,ierror,jacobn,matdim,iwork(imxord),iwork(imnt),iwork(imtr),ml,mu,npar,ndecom,work(sub2ind(size(work),max(iywt,1)):end),uround,users,work(iavgh),work(iavgrd),work(ih),hused,iwork(ijtask),iwork(imntld),iwork(imtrld),iwork(infe),iwork(inje),iwork(inquse),iwork(instep),work(it),y,work(sub2ind(size(work),max(iyh,1)):end),work(sub2ind(size(work),max(ia,1)):end),convrg,work(sub2ind(size(work),max(idfdy,1)):end),work(iel:end),work(sub2ind(size(work),max(ifac,1)):end),work(ihold),iwork(sub2ind(size(iwork),max(indpvt,1)):end),jstate,iwork(ijstpl),iwork(inq),iwork(inwait),work(irc),work(irmax),work(sub2ind(size(work),max(isave1,1)):end),work(sub2ind(size(work),max(isave2,1)):end),work(itq:end),work(itrend),mint,iwork(imtrsv),iwork(imxrds)); dumvar4i=find((work(ihmax))~=(dumvar4));dumvar9i=find((iwork(imxord))~=(dumvar9));dumvar10i=find((iwork(imnt))~=(dumvar10));dumvar11i=find((iwork(imtr))~=(dumvar11));dumvar16i=find((work(sub2ind(size(work),max(iywt,1)):end))~=(dumvar16));dumvar19i=find((work(iavgh))~=(dumvar19));dumvar20i=find((work(iavgrd))~=(dumvar20));dumvar21i=find((work(ih))~=(dumvar21));dumvar23i=find((iwork(ijtask))~=(dumvar23));dumvar24i=find((iwork(imntld))~=(dumvar24));dumvar25i=find((iwork(imtrld))~=(dumvar25));dumvar26i=find((iwork(infe))~=(dumvar26));dumvar27i=find((iwork(inje))~=(dumvar27));dumvar28i=find((iwork(inquse))~=(dumvar28));dumvar29i=find((iwork(instep))~=(dumvar29));dumvar30i=find((work(it))~=(dumvar30));dumvar32i=find((work(sub2ind(size(work),max(iyh,1)):end))~=(dumvar32));dumvar33i=find((work(sub2ind(size(work),max(ia,1)):end))~=(dumvar33));dumvar35i=find((work(sub2ind(size(work),max(idfdy,1)):end))~=(dumvar35));dumvar36i=find((work(iel:end))~=(dumvar36));dumvar37i=find((work(sub2ind(size(work),max(ifac,1)):end))~=(dumvar37));dumvar38i=find((work(ihold))~=(dumvar38));dumvar39i=find((iwork(sub2ind(size(iwork),max(indpvt,1)):end))~=(dumvar39));dumvar41i=find((iwork(ijstpl))~=(dumvar41));dumvar42i=find((iwork(inq))~=(dumvar42));dumvar43i=find((iwork(inwait))~=(dumvar43));dumvar44i=find((work(irc))~=(dumvar44));dumvar45i=find((work(irmax))~=(dumvar45));dumvar46i=find((work(sub2ind(size(work),max(isave1,1)):end))~=(dumvar46));dumvar47i=find((work(sub2ind(size(work),max(isave2,1)):end))~=(dumvar47));dumvar48i=find((work(itq:end))~=(dumvar48));dumvar49i=find((work(itrend))~=(dumvar49));dumvar51i=find((iwork(imtrsv))~=(dumvar51));dumvar52i=find((iwork(imxrds))~=(dumvar52)); work(ihmax-1+dumvar4i)=dumvar4(dumvar4i); iwork(imxord-1+dumvar9i)=dumvar9(dumvar9i); iwork(imnt-1+dumvar10i)=dumvar10(dumvar10i); iwork(imtr-1+dumvar11i)=dumvar11(dumvar11i); work(iywt-1+dumvar16i)=dumvar16(dumvar16i); work(iavgh-1+dumvar19i)=dumvar19(dumvar19i); work(iavgrd-1+dumvar20i)=dumvar20(dumvar20i); work(ih-1+dumvar21i)=dumvar21(dumvar21i); iwork(ijtask-1+dumvar23i)=dumvar23(dumvar23i); iwork(imntld-1+dumvar24i)=dumvar24(dumvar24i); iwork(imtrld-1+dumvar25i)=dumvar25(dumvar25i); iwork(infe-1+dumvar26i)=dumvar26(dumvar26i); iwork(inje-1+dumvar27i)=dumvar27(dumvar27i); iwork(inquse-1+dumvar28i)=dumvar28(dumvar28i); iwork(instep-1+dumvar29i)=dumvar29(dumvar29i); work(it-1+dumvar30i)=dumvar30(dumvar30i); work(iyh-1+dumvar32i)=dumvar32(dumvar32i); work(ia-1+dumvar33i)=dumvar33(dumvar33i); work(idfdy-1+dumvar35i)=dumvar35(dumvar35i); work(iel-1+dumvar36i)=dumvar36(dumvar36i); work(ifac-1+dumvar37i)=dumvar37(dumvar37i); work(ihold-1+dumvar38i)=dumvar38(dumvar38i); iwork(indpvt-1+dumvar39i)=dumvar39(dumvar39i); iwork(ijstpl-1+dumvar41i)=dumvar41(dumvar41i); iwork(inq-1+dumvar42i)=dumvar42(dumvar42i); iwork(inwait-1+dumvar43i)=dumvar43(dumvar43i); work(irc-1+dumvar44i)=dumvar44(dumvar44i); work(irmax-1+dumvar45i)=dumvar45(dumvar45i); work(isave1-1+dumvar46i)=dumvar46(dumvar46i); work(isave2-1+dumvar47i)=dumvar47(dumvar47i); work(itq-1+dumvar48i)=dumvar48(dumvar48i); work(itrend-1+dumvar49i)=dumvar49(dumvar49i); iwork(imtrsv-1+dumvar51i)=dumvar51(dumvar51i); iwork(imxrds-1+dumvar52i)=dumvar52(dumvar52i);
t = work(it);
h = work(ih);
if( convrg )
iwork(icnvrg) = 1;
else;
iwork(icnvrg) = 0;
end;
if( jstate==2 )
gt(6)=1;
continue;
end;
if( jstate==3 )
%
rl1=sprintf(['%16.8f'], t);
ierflg = 42;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,[', the step sizemlv has been reduced about 50 ',['times without advancing the solution. Often this occurs ','if the problem setup is incorrect.']]]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
elseif( jstate==4 || jstate==5 ) ;
break;
elseif( jstate==6 || jstate==7 || jstate==8 ||jstate==9 || jstate==10 ) ;
%
nstate = fix(jstate);
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
else;
iwork(ijtask) = 1;
% Determine if a root has been overtaken
if( nroot~=0 )
iroot = 0;
for i = 1 : nroot;
glast = work(i+ignow-1);
gnow = g(npar,t,y,i);
if( npar==0 )
iwork(inroot) = fix(i);
nstate = 7;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
work(i+ignow-1) = gnow;
if( glast.*gnow>0.0e0 )
work(i+itroot-1) = t + h;
elseif( gnow==0.0e0 ) ;
work(i+itroot-1) = t;
iroot = fix(i);
elseif( glast==0.0e0 ) ;
work(i+itroot-1) = t + h;
elseif( abs(hused)>=uround.*abs(t) ) ;
tlast = t - hused;
iroot = fix(i);
troot = t;
[ae,g,h,npar,iwork(inq),iroot,re,t,work(sub2ind(size(work),max(iyh,1)):end),uround,troot,tlast,gnow,glast,y]=sdzro(ae,g,h,npar,iwork(inq),iroot,re,t,work(sub2ind(size(work),max(iyh,1)):end),uround,troot,tlast,gnow,glast,y);
for j = 1 : n;
y(j) = work(iyh+j-1);
end; j = fix(n+1);
if( npar==0 )
iwork(inroot) = fix(i);
nstate = 7;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
work(i+itroot-1) = troot;
else;
work(i+itroot-1) = t;
iroot = fix(i);
end;
end; i = fix(nroot+1);
if( iroot==0 )
iwork(ijroot) = 0;
% Select the first root
else;
iwork(ijroot) = fix(ntask);
iwork(inrtld) = fix(nroot);
iwork(indtrt) = fix(itroot);
troot = t + h;
for i = 1 : nroot;
if( work(i+itroot-1).*hsignmlv<troot.*hsignmlv )
troot = work(i+itroot-1);
iroot = fix(i);
end;
end; i = fix(nroot+1);
iwork(inroot) = fix(iroot);
work(itout) = troot;
if( troot.*hsignmlv<=tout.*hsignmlv )
[h,dumvar2,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,troot,work(sub2ind(size(work),max(iyh,1)):end),y);
nstate = 5;
t = troot;
ierflg = 0;
gt(5)=1;
continue;
end;
end;
end;
% Test for NTASK condition to be satisfied
nstate = 2;
if( ntask==1 )
if( t.*hsignmlv<tout.*hsignmlv )
gt(2)=1;
continue;
end;
[h,dumvar2,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y]=sdntp(h,0,n,iwork(inq),t,tout,work(sub2ind(size(work),max(iyh,1)):end),y);
t = tout;
ierflg = 0;
gt(5)=1;
continue;
% TOUT is assumed to have been attained
% exactly if T is within twenty roundoff
% units of TOUT, relative to MAX(TOUT, T).
%
elseif( ntask==2 ) ;
if( abs(tout-t)<=nround.*uround.*max(abs(t),abs(tout)) )
t = tout;
elseif((t+h).*hsignmlv>tout.*hsignmlv ) ;
h = tout - t;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = h.*(1.0e0-4.0e0.*uround);
end;
work(ih) = h;
if( h==0.0e0 )
gt(6)=1;
continue;
end;
iwork(ijtask) = -1;
end;
elseif( ntask==3 ) ;
if( abs(tout-t)<=nround.*uround.*max(abs(t),abs(tout)) )
t = tout;
else;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = tout - t;
if((t+h).*hsignmlv>tout.*hsignmlv )
h = h.*(1.0e0-4.0e0.*uround);
end;
work(ih) = h;
if( h==0.0e0 )
gt(6)=1;
continue;
end;
iwork(ijtask) = -1;
end;
gt(2)=1;
continue;
end;
end;
ierflg = 0;
end;
% All returns are made through this
% section. IMXERR is determined.
end;
gt(4)=0;
for i = 1 : n;
y(i) = work(i+iyh-1);
end; i = fix(n+1);
end;
gt(5)=0;
if( iwork(ijtask)==0 )
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
end;
big = 0.0e0;
imxerr = 1;
for i = 1 : n;
% SIZE = ABS(ERROR(I)/YWT(I))
sizemlv = abs(work(i+isave1-1)./work(i+iywt-1));
if( big<sizemlv )
big = sizemlv;
imxerr = fix(i);
end;
end; i = fix(n+1);
iwork(indmxr) = fix(imxerr);
work(ihused) = hused;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
% Fatal errors are processed here
%
end;
gt(6)=0;
rl1=sprintf(['%16.8f'], t);
ierflg = 41;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,[', the attempted step sizemlv has gone to ','zero. Often this occurs if the problem setup is incorrect.']]],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
return;
%
break;
end;
rl1=sprintf(['%16.8f'], t);
ierflg = 43;
[dumvar1,dumvar2,dumvar3,ierflg]=xermsg('SLATEC','SDRIV3',['At T, ',[rl1,', while solving A*YDOT = F, A is singular.']],ierflg,1);
nstate = 12;
ewt_shape=zeros(ewt_shape);ewt_shape(:)=ewt(1:numel(ewt_shape));ewt=ewt_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
end %subroutine sdriv3
%DECK SDSCL
|
|
