Math-Lsoda
view release on metacpan or search on metacpan
C
C Note: The work arrays must not be altered between calls to DLSODE
C for the same problem, except possibly for the conditional and
C optional inputs, and except for the last 3*NEQ words of RWORK.
C The latter space is used for internal scratch space, and so is
C available for use by the user outside DLSODE between calls, if
C desired (but not for use by F or JAC).
C
C JAC The name of the user-supplied routine (MITER = 1 or 4) to
C compute the Jacobian matrix, df/dy, as a function of the
C scalar t and the vector y. (See the MF description below
C for MITER.) It is to have the form
C
C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
C DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
C
C where NEQ, T, Y, ML, MU, and NROWPD are input and the
C array PD is to be loaded with partial derivatives
C (elements of the Jacobian matrix) on output. PD must be
C given a first dimension of NROWPD. T and Y have the same
C meaning as in subroutine F.
C
C In the full matrix case (MITER = 1), ML and MU are
C ignored, and the Jacobian is to be loaded into PD in
C columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
C
C In the band matrix case (MITER = 4), the elements within
C the band are to be loaded into PD in columnwise manner,
C with diagonal lines of df/dy loaded into the rows of PD.
C Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML
C and MU are the half-bandwidth parameters (see IWORK).
C The locations in PD in the two triangular areas which
C correspond to nonexistent matrix elements can be ignored
C or loaded arbitrarily, as they are overwritten by DLSODE.
C
C JAC need not provide df/dy exactly. A crude approximation
C (possibly with a smaller bandwidth) will do.
C
C In either case, PD is preset to zero by the solver, so
C that only the nonzero elements need be loaded by JAC.
C Each call to JAC is preceded by a call to F with the same
C arguments NEQ, T, and Y. Thus to gain some efficiency,
C intermediate quantities shared by both calculations may
C be saved in a user COMMON block by F and not recomputed
C by JAC, if desired. Also, JAC may alter the Y array, if
C desired. JAC must be declared EXTERNAL in the calling
C program.
C
C Subroutine JAC may access user-defined quantities in
C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in JAC) and/or Y has length exceeding
C NEQ(1). See the descriptions of NEQ and Y above.
C
C MF The method flag. Used only for input. The legal values
C of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24,
C and 25. MF has decimal digits METH and MITER:
C MF = 10*METH + MITER .
C
C METH indicates the basic linear multistep method:
C 1 Implicit Adams method.
C 2 Method based on backward differentiation formulas
C (BDF's).
C
C MITER indicates the corrector iteration method:
C 0 Functional iteration (no Jacobian matrix is
C involved).
C 1 Chord iteration with a user-supplied full (NEQ by
C NEQ) Jacobian.
C 2 Chord iteration with an internally generated
C (difference quotient) full Jacobian (using NEQ
C extra calls to F per df/dy value).
C 3 Chord iteration with an internally generated
C diagonal Jacobian approximation (using one extra call
C to F per df/dy evaluation).
C 4 Chord iteration with a user-supplied banded Jacobian.
C 5 Chord iteration with an internally generated banded
C Jacobian (using ML + MU + 1 extra calls to F per
C df/dy evaluation).
C
C If MITER = 1 or 4, the user must supply a subroutine JAC
C (the name is arbitrary) as described above under JAC.
C For other values of MITER, a dummy argument can be used.
C
C Optional Inputs
C ---------------
C The following is a list of the optional inputs provided for in the
C call sequence. (See also Part 2.) For each such input variable,
C this table lists its name as used in this documentation, its
C location in the call sequence, its meaning, and the default value.
C The use of any of these inputs requires IOPT = 1, and in that case
C all of these inputs are examined. A value of zero for any of
C these optional inputs will cause the default value to be used.
C Thus to use a subset of the optional inputs, simply preload
C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively,
C and then set those of interest to nonzero values.
C
C Name Location Meaning and default value
C ------ --------- -----------------------------------------------
C H0 RWORK(5) Step size to be attempted on the first step.
C The default value is determined by the solver.
C HMAX RWORK(6) Maximum absolute step size allowed. The
C default value is infinite.
C HMIN RWORK(7) Minimum absolute step size allowed. The
C default value is 0. (This lower bound is not
C enforced on the final step before reaching
C TCRIT when ITASK = 4 or 5.)
C MAXORD IWORK(5) Maximum order to be allowed. The default value
C is 12 if METH = 1, and 5 if METH = 2. (See the
C MF description above for METH.) If MAXORD
C exceeds the default value, it will be reduced
C to the default value. If MAXORD is changed
C during the problem, it may cause the current
C order to be reduced.
C MXSTEP IWORK(6) Maximum number of (internally defined) steps
C allowed during one call to the solver. The
C default value is 500.
C MXHNIL IWORK(7) Maximum number of messages printed (per
C problem) warning that T + H = T on a step
C (H = step size). This must be positive to
C result in a nondefault value. The default
C value is 10.
C to 3. The optional output TOLSF may be used for this
C purpose. (Note: If this condition is detected before
C taking any steps, then an illegal input return
C (ISTATE = -3) occurs instead.)
C -3 means illegal input was detected, before taking any
C integration steps. See written message for details.
C Note: If the solver detects an infinite loop of calls
C to the solver with illegal input, it will cause
C the run to stop.
C -4 means there were repeated error test failures on
C one attempted step, before completing the requested
C task, but the integration was successful as far as T.
C The problem may have a singularity, or the input
C may be inappropriate.
C -5 means there were repeated convergence test failures on
C one attempted step, before completing the requested
C task, but the integration was successful as far as T.
C This may be caused by an inaccurate Jacobian matrix,
C if one is being used.
C -6 means EWT(i) became zero for some i during the
C integration. Pure relative error control (ATOL(i)=0.0)
C was requested on a variable which has now vanished.
C The integration was successful as far as T.
C -7 means a fatal error return flag came from the sparse
C solver CDRV by way of DPRJS or DSOLSS (numerical
C factorization or backsolve). This should never happen.
C The integration was successful as far as T.
C
C Note: an error return with ISTATE = -1, -4, or -5 and with
C MITER = 1 or 2 may mean that the sparsity structure of the
C problem has changed significantly since it was last
C determined (or input). In that case, one can attempt to
C complete the integration by setting ISTATE = 3 on the next
C call, so that a new structure determination is done.
C
C Note: since the normal output value of ISTATE is 2,
C it does not need to be reset for normal continuation.
C Also, since a negative input value of ISTATE will be
C regarded as illegal, a negative output value requires the
C user to change it, and possibly other inputs, before
C calling the solver again.
C
C IOPT = an integer flag to specify whether or not any optional
C inputs are being used on this call. Input only.
C The optional inputs are listed separately below.
C IOPT = 0 means no optional inputs are being used.
C Default values will be used in all cases.
C IOPT = 1 means one or more optional inputs are being used.
C
C RWORK = a work array used for a mixture of real (double precision)
C and integer work space.
C The length of RWORK (in real words) must be at least
C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM where
C NYH = the initial value of NEQ,
C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
C smaller value is given as an optional input),
C LWM = 0 if MITER = 0,
C LWM = 2*NNZ + 2*NEQ + (NNZ+9*NEQ)/LENRAT if MITER = 1,
C LWM = 2*NNZ + 2*NEQ + (NNZ+10*NEQ)/LENRAT if MITER = 2,
C LWM = NEQ + 2 if MITER = 3.
C In the above formulas,
C NNZ = number of nonzero elements in the Jacobian matrix.
C LENRAT = the real to integer wordlength ratio (usually 1 in
C single precision and 2 in double precision).
C (See the MF description for METH and MITER.)
C Thus if MAXORD has its default value and NEQ is constant,
C the minimum length of RWORK is:
C 20 + 16*NEQ for MF = 10,
C 20 + 16*NEQ + LWM for MF = 11, 111, 211, 12, 112, 212,
C 22 + 17*NEQ for MF = 13,
C 20 + 9*NEQ for MF = 20,
C 20 + 9*NEQ + LWM for MF = 21, 121, 221, 22, 122, 222,
C 22 + 10*NEQ for MF = 23.
C If MITER = 1 or 2, the above formula for LWM is only a
C crude lower bound. The required length of RWORK cannot
C be readily predicted in general, as it depends on the
C sparsity structure of the problem. Some experimentation
C may be necessary.
C
C The first 20 words of RWORK are reserved for conditional
C and optional inputs and optional outputs.
C
C The following word in RWORK is a conditional input:
C RWORK(1) = TCRIT = critical value of t which the solver
C is not to overshoot. Required if ITASK is
C 4 or 5, and ignored otherwise. (See ITASK.)
C
C LRW = the length of the array RWORK, as declared by the user.
C (This will be checked by the solver.)
C
C IWORK = an integer work array. The length of IWORK must be at least
C 31 + NEQ + NNZ if MOSS = 0 and MITER = 1 or 2, or
C 30 otherwise.
C (NNZ is the number of nonzero elements in df/dy.)
C
C In DLSODES, IWORK is used only for conditional and
C optional inputs and optional outputs.
C
C The following two blocks of words in IWORK are conditional
C inputs, required if MOSS = 0 and MITER = 1 or 2, but not
C otherwise (see the description of MF for MOSS).
C IWORK(30+j) = IA(j) (j=1,...,NEQ+1)
C IWORK(31+NEQ+k) = JA(k) (k=1,...,NNZ)
C The two arrays IA and JA describe the sparsity structure
C to be assumed for the Jacobian matrix. JA contains the row
C indices where nonzero elements occur, reading in columnwise
C order, and IA contains the starting locations in JA of the
C descriptions of columns 1,...,NEQ, in that order, with
C IA(1) = 1. Thus, for each column index j = 1,...,NEQ, the
C values of the row index i in column j where a nonzero
C element may occur are given by
C i = JA(k), where IA(j) .le. k .lt. IA(j+1).
C If NNZ is the total number of nonzero locations assumed,
C then the length of the JA array is NNZ, and IA(NEQ+1) must
C be NNZ + 1. Duplicate entries are not allowed.
C
C LIW = the length of the array IWORK, as declared by the user.
C (This will be checked by the solver.)
C
C Note: The work arrays must not be altered between calls to DLSODES
C for the same problem, except possibly for the conditional and
C optional inputs, and except for the last 3*NEQ words of RWORK.
C The latter space is used for internal scratch space, and so is
C available for use by the user outside DLSODES between calls, if
C desired (but not for use by F or JAC).
C
C JAC = name of user-supplied routine (MITER = 1 or MOSS = 1) to
C compute the Jacobian matrix, df/dy, as a function of
C the scalar t and the vector y. It is to have the form
C SUBROUTINE JAC (NEQ, T, Y, J, IAN, JAN, PDJ)
C DOUBLE PRECISION T, Y(*), IAN(*), JAN(*), PDJ(*)
C where NEQ, T, Y, J, IAN, and JAN are input, and the array
C PDJ, of length NEQ, is to be loaded with column J
C of the Jacobian on output. Thus df(i)/dy(J) is to be
C Note: If the solver detects an infinite loop of calls
C to the solver with illegal input, it will cause
C the run to stop.
C -4 means there were repeated error test failures on
C one attempted step, before completing the requested
C task, but the integration was successful as far as T.
C The problem may have a singularity, or the input
C may be inappropriate.
C -5 means there were repeated convergence test failures on
C one attempted step, before completing the requested
C task, but the integration was successful as far as T.
C This may be caused by an inaccurate Jacobian matrix.
C -6 means EWT(i) became zero for some i during the
C integration. Pure relative error control (ATOL(i) = 0.0)
C was requested on a variable which has now vanished.
C the integration was successful as far as T.
C -7 means that the user-supplied Subroutine RES set
C its error flag (IRES = 3) despite repeated tries by
C DLSODIS to avoid that condition.
C -8 means that ISTATE was 0 on input but DLSODIS was unable
C to compute the initial value of dy/dt. See the
C printed message for details.
C -9 means a fatal error return flag came from the sparse
C solver CDRV by way of DPRJIS or DSOLSS (numerical
C factorization or backsolve). This should never happen.
C The integration was successful as far as T.
C
C Note: An error return with ISTATE = -1, -4, or -5
C may mean that the sparsity structure of the
C problem has changed significantly since it was last
C determined (or input). In that case, one can attempt to
C complete the integration by setting ISTATE = 3 on the next
C call, so that a new structure determination is done.
C
C Note: Since the normal output value of ISTATE is 2,
C it does not need to be reset for normal continuation.
C similarly, ISTATE (= 3) need not be reset if RES told
C DLSODIS to return because the calling program must change
C the parameters of the problem.
C Also, since a negative input value of ISTATE will be
C regarded as illegal, a negative output value requires the
C user to change it, and possibly other inputs, before
C calling the solver again.
C
C IOPT = an integer flag to specify whether or not any optional
C inputs are being used on this call. Input only.
C The optional inputs are listed separately below.
C IOPT = 0 means no optional inputs are being used.
C Default values will be used in all cases.
C IOPT = 1 means one or more optional inputs are being used.
C
C RWORK = a work array used for a mixture of real (double precision)
C and integer work space.
C The length of RWORK (in real words) must be at least
C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM where
C NYH = the initial value of NEQ,
C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
C smaller value is given as an optional input),
C LWM = 2*NNZ + 2*NEQ + (NNZ+9*NEQ)/LENRAT if MITER = 1,
C LWM = 2*NNZ + 2*NEQ + (NNZ+10*NEQ)/LENRAT if MITER = 2.
C in the above formulas,
C NNZ = number of nonzero elements in the iteration matrix
C P = A - con*J (con is a constant and J is the
C Jacobian matrix dr/dy).
C LENRAT = the real to integer wordlength ratio (usually 1 in
C single precision and 2 in double precision).
C (See the MF description for METH and MITER.)
C Thus if MAXORD has its default value and NEQ is constant,
C the minimum length of RWORK is:
C 20 + 16*NEQ + LWM for MF = 11, 111, 311, 12, 212, 412,
C 20 + 9*NEQ + LWM for MF = 21, 121, 321, 22, 222, 422.
C The above formula for LWM is only a crude lower bound.
C The required length of RWORK cannot be readily predicted
C in general, as it depends on the sparsity structure
C of the problem. Some experimentation may be necessary.
C
C The first 20 words of RWORK are reserved for conditional
C and optional inputs and optional outputs.
C
C The following word in RWORK is a conditional input:
C RWORK(1) = TCRIT = critical value of t which the solver
C is not to overshoot. Required if ITASK is
C 4 or 5, and ignored otherwise. (See ITASK.)
C
C LRW = the length of the array RWORK, as declared by the user.
C (This will be checked by the solver.)
C
C IWORK = an integer work array. The length of IWORK must be at least
C 32 + 2*NEQ + NZA + NZC for MOSS = 0,
C 30 for MOSS = 1 or 2,
C 31 + NEQ + NZA for MOSS = 3 or 4.
C (NZA is the number of nonzero elements in matrix A, and
C NZC is the number of nonzero elements in dr/dy.)
C
C In DLSODIS, IWORK is used for conditional and
C optional inputs and optional outputs.
C
C The following two blocks of words in IWORK are conditional
C inputs, required if MOSS = 0, 3, or 4, but not otherwise
C (see the description of MF for MOSS).
C IWORK(30+j) = IA(j) (j=1,...,NEQ+1)
C IWORK(31+NEQ+k) = JA(k) (k=1,...,NZA)
C The two arrays IA and JA describe the sparsity structure
C to be assumed for the matrix A. JA contains the row
C indices where nonzero elements occur, reading in columnwise
C order, and IA contains the starting locations in JA of the
C descriptions of columns 1,...,NEQ, in that order, with
C IA(1) = 1. Thus, for each column index j = 1,...,NEQ, the
C values of the row index i in column j where a nonzero
C element may occur are given by
C i = JA(k), where IA(j) .le. k .lt. IA(j+1).
C If NZA is the total number of nonzero locations assumed,
C then the length of the JA array is NZA, and IA(NEQ+1) must
C be NZA + 1. Duplicate entries are not allowed.
C The following additional blocks of words are required
C if MOSS = 0, but not otherwise. If LC = 31 + NEQ + NZA, then
C IWORK(LC+j) = IC(j) (j=1,...,NEQ+1), and
C IWORK(LC+NEQ+1+k) = JC(k) (k=1,...,NZC)
C The two arrays IC and JC describe the sparsity
C structure to be assumed for the Jacobian matrix dr/dy.
C They are used in the same manner as the above IA and JA
C arrays. If NZC is the number of nonzero locations
C assumed, then the length of the JC array is NZC, and
C IC(NEQ+1) must be NZC + 1. Duplicate entries are not
C allowed.
C
C LIW = the length of the array IWORK, as declared by the user.
C (This will be checked by the solver.)
C
C Note: The work arrays must not be altered between calls to DLSODIS
C for the same problem, except possibly for the conditional and
C optional inputs, and except for the last 3*NEQ words of RWORK.
( run in 0.534 second using v1.01-cache-2.11-cpan-39bf76dae61 )