DAESL

more...

Solves a first order differential-algebraic system of equations, g(t, y, yʹ) = 0, with optional additional constraints and user-defined linear system solver.

Note: DAESL replaces deprecated routine DASPG.

Required Arguments

T — Independent variable, t. (Input/Output)
Set T to the starting value t0 at the first step. On output, T is set to the value to which the integration has advanced. Normally, this new value is TEND.

TEND — Final value of the independent variable. (Input)
Update this value when re-entering after output with IDO = 2.

IDO — Flag indicating the state of the computation. (Input/Output)

IDO	State
1	Initial entry
2	Normal re-entry after obtaining output
3	Release workspace, last call

The user sets IDO = 1 on the first call at T = t0. The routine then sets IDO = 2, and this value is used for all but the last entry, which is made with IDO = 3.

Y — Array of size NEQ containing the dependent variable values, y. (Input/Output)
On input, Y must contain initial values. On output, Y contains the computed solution at TEND.

YPRIME — Array of size NEQ containing derivative values, yʹ. (Input/Output)
This array must contain initial values, but they need not be such that g(t, y, yʹ) = 0 at t = t0. See the description of parameter IYPR for more information.

GCN — User-supplied subroutine to evaluate g(t, y, yʹ), and any constraints. Also partial derivative evaluations and optionally linear solving steps occur here. The equations g(t, y, yʹ) = 0 consist of NEQ differential-algebraic equations of the form.

Fi(t,y1, … yNEQ, y1′, … y′NEQ) ≡ Fi(t,y,y′) = 0, i = 1, …, NEQ

The routine GCN is also used to evaluate the NCON additional algebraic constraints

Gi(t, y1, …, yNEQ) ≡ Gi(t, y) = 0, i = 1, …, NCONNCON ≥ 0

The usage is CALL GCN (T, Y, YPRIME, DELTA, D, LDD, IRES [,…]) where

Required Arguments

T — Integration variable t. (Input)

Y — Array of NEQ dependent variables, y. (Input)

YPRIME — Array of NEQ derivative values, yʹ. (Input)

DELTA — Output array of length MAX(NEQ, NCON) containing residuals. See parameter IRES for definition. (Input/Output)

D — Output array dimensioned D(LDD,NEQ), containing partial derivatives. See parameter IRES for definition. (Input/Output)

LDD — Leading dimension of D. (Input)

IRES — Flag indicating what is to be calculated in the user routine, GCN. (Input/Output)
Note: IRES is input only, except when IRES = 6. It is input/output when IRES = 6. For a detailed description see the table below.

The code calls GCN with IRES = 0, 1, 2, 3, 4, 5, 6, or 7, defined as follows:

IRES Value	Explanation
0	Do initializations, if any are required.
1	Compute DELTA(i) = , the i-th residual, for i = 1, …, NEQ.
2	(Required only if IUJAC = 1 and MATSTR = 0 or 1.) Compute D(i, j) = , the partial derivative matrix. These are derivatives of Fi with respect to yj, for i = 1, …, NEQ and j = 1, …, NEQ.
3	(Required only if IUJAC = 1 and MATSTR = 0 or 1.) Compute D(i, j) = , the partial derivative matrix. These are derivatives of Fi with respect to yj′, for i = 1, …, NEQ and j = 1, …, NEQ.
4	(Required only if IYPR=2.) Compute DELTA(i) = , the partial derivative of Fi with respect to t, for i = 1, …, NEQ.
5	(Required only if NCON > 0.) Compute DELTA(i) = Gi(t, y), the i-th residual in the additional constraints, for i = 1, …, NCON, and D(i, j) = , the partial derivative of Gi with respect to yj for i = 1, …, NCON and j = 1, …, NEQ.
6	(Required only if ISOLVE = 1.) If MATSTR = 2, the user must compute the matrix , where cj = DELTA (1), and save this matrix in any user-defined format. This is for later use when IRES = 7. The matrix may also be factored in this step, if desired. The array D is not referenced if MATSTR = 2. If MATSTR = 0 or 1, the A matrix will already be defined and passed to GCN in the array D, which will be in full matrix format if MATSTR = 0, and band matrix format, if MATSTR = 1. The user may factor D in this step, if desired. Note: For MATSTR = 0, 1, or 2, the user must set IRES = 0 to signal that A is nonsingular. If A is nearly singular, leave IRES = 6. This results in using a smaller step-size internally.
7	(Required only if ISOLVE = 1.) The user must solve Ax = b, where b is passed to GCN in the vector DELTA, and x is returned in DELTA. If MATSTR = 2, A is the matrix which was computed and saved at the call with IRES = 6; if MATSTR = 0 or 1, A is passed to GCN in the array D. In either case, the A matrix will remain factored if the user factored it when IRES = 6.

Optional Arguments

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a description of this argument see FCN_DATA below. (Input/Output)

GCN must be declared EXTERNAL in the calling program.

Optional Arguments

NEQ — Number of dependent variables, and number of differential/algebraic equations, not counting any additional constraints. (Input)
Default: NEQ = size (Y).

NCON — Number of additional constraints. (Input)
Default: NCON = 0.

IUJAC — Jacobian calculation option. (Input)

Value	Description
0	Calculates using finite difference approximations.
1	User supplies the Jacobian matrices of partial derivatives of Fi, i = 1, …, NEQ, in the subroutine GCN, when IRES = 2 and 3.

Default: IUJAC = 0 for MATSTR = 0 or 1.
IUJAC = 1 for MATSTR = 2.

IYPR — Initial yʹ calculation method. (Input)

Value	Description
0	The initial input values of YPRIME are already consistent with the input values of Y. That is g(t, y, yʹ) = 0 at t = t0. Any constraints must be satisfied at t = t0.
1	Consistent values of YPRIME are calculated by Petzold’s original DASSL algorithm.
2	Consistent values of YPRIME are calculated using a new algorithm [Hanson and Krogh, 2008], which is generally more robust but requires that IUJAC = 1 and ISOLVE = 0, and additional derivatives corresponding to IRES = 4 are to be calculated in GCN.

Default: IYPR = 1.

MATSTR — Parameter specifying the Jacobian matrix structure (Input)
Set to:

Value	Description
0	The Jacobian matrices (whether IUJAC = 0 or 1) are to be stored in full storage mode.
1	The Jacobian matrices are to be stored in band storage mode. In this case, if IUJAC = 1, the partial derivative matrices have their entries for row i and column j, stored as array elements D(i ‑ j + MU + 1, j). This occurs when IRES= 2 or 3 in GCN.
2	A user-defined matrix structure is used (see the documentation for IRES = 6 or 7 for more details). If MATSTR = 2, ISOLVE and IUJAC are set to 1 internally.

Default: MATSTR = 0.

ISOLVE — Solve method. (Input)

Value	Description
0	DAESL solves the linear systems.
1	The user wishes to solve the linear system in routine GCN. See parameter GCN for details.

Default: ISOLVE = 0 for MATSTR = 0 or 1, ISOLVE = 1 for MATSTR = 2.

ML — Number of non-zero diagonals below the main diagonal in the Jacobian matrices when band storage mode is used. (Input)
ML is ignored if MATSTR ≠ 1.
Default: ML = NEQ-1.

MU — Number of non-zero diagonals above the main diagonal in the Jacobian matrices when band storage mode is used. (Input)
MU is ignored if MATSTR ≠ 1.
Default: MU = NEQ-1.

RTOL — Relative error tolerance for solver. (Input)
The program attempts to maintain a local error in Y(i) less than
RTOL*∣Y(i)∣ + ATOL(i).
Default: RTOL = , where ɛ is machine precision.

ATOL — Array of size NEQ containing absolute error tolerances. (Input)
See description of RTOL.
Default: ATOL(i) = 0.

H0 — Initial stepsize used by the solver. (Input)
If H0 = 0, the routine defines the initial stepsize.
Default: H0 = 0.

HMAX — Maximum stepsize used by the solver. (Input)
If HMAX=0, the routine defines the maximum stepsize.
Default: HMAX = 0.

MAXORD — Maximum order of the backward difference formulas used. (Input).
1 ≤ MAXORD ≤ 5.
Default: MAXORD = 5.

MAXSTEPS — Maximum number of steps taken from T to TEND. (Input).
Default: MAXSTEPS = 500.

TSTOP — Integration limit point. (Input)
For efficiency reasons, the code sometimes integrates past TEND and interpolates a solution at TEND. If a value for TSTOP is specified, the code will never integrate past T=TSTOP.
Default: No TSTOP value is specified.

FMAG — Order-of-magnitude estimate. (Input)
FMAG is used as an order-of-magnitude estimate of the magnitude of the functions Fi (see description of GCN), for convergence testing, if IYPR = 2. FMAG is ignored if IYPR=0 or 1.
Default: FMAG = 1.

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied subroutine. (Input/Output)
The derived type, s_fcn_data, is defined as:

type s_fcn_data

real(kind(1e0)), pointer, dimension(:) :: rdata

integer, pointer, dimension(:) :: idata

end type

in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type.

Note that if this optional argument is present then FCN_DATA must also be defined as an optional argument in the user-supplied subroutine.

FORTRAN 90 Interface

Generic: CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN [, …])

Specific: The specific interface names are S_DAESL and D_DAESL.

Description

Routine DAESL finds an approximation to the solution of a system of differential-algebraic equations g(t, y, y′) = 0 with given initial data for y and y′. The routine uses BDF formulas, which are appropriate for stiff systems. DAESL is based on the code DASSL designed by Linda Petzold [1982], and has been modified by Hanson and Krogh [2008] Solving Constrained Differential-Algebraic Systems Using Projections to allow the inclusion of additional constraints, including conservation principles, after each time step. The modified code also provides a more robust algorithm to calculate initial y′ values consistent with the given initial y values. This occurs when the initial y′ are not known.

A differential-algebraic system of equations is said to have “index 0” if the Jacobian matrix of partial derivatives of the Fi with respect to the yj′ is nonsingular. Thus it is possible to solve for all the initial values of yj′ and put the system in the form of a standard ODE system. If it is possible to reduce the system to a system of index 0 by taking first derivatives of some of the equations, the system has index 1, otherwise the index is greater than 1. See Brenan [1989] for a definition of index. DAESL can generally only solve systems of index 0 or 1; other systems will usually have to be reduced to such a form through differentiation.

Examples

Example 1 – Method of Lines PDE Problem

This example solves the partial differential equation , with initial condition , and boundary conditions , which has exact solution . If we approximate the term using finite differences, where

, and , we get:

If Yi(t) = U(xi,t), the first and last equations are algebraic and the others are differential equations, so this is a system of differential-algebraic equations. The system has index = 1, since it could be transformed into an ODE system by differentiating the first and last equations. Note that the Jacobian matrices are banded (tridiagonal), with ML = MU = 1. We use this and specify the option for dealing with banded matrices in DAESL. The parameter and the number of equations is passed to the evaluation routine, GCN, with the optional argument USER_DATA.

(Example daesl_ex1.f90)

USE DAESL_INT

USE MP_TYPES

IMPLICIT NONE

! NEQ = Number of equations

INTEGER, PARAMETER :: NEQ=101

REAL T, Y(NEQ), YPRIME(NEQ), TEND, X, TRUE, HX, ERRMAX

INTEGER NOUT, IDO, I, NSTEPS

REAL, TARGET :: RPARAM(1)

INTEGER, TARGET :: IPARAM(1)

TYPE (S_FCN_DATA) USER_DATA

EXTERNAL GCN

! Pass NEQ, HX to GCN

HX = 1.0 / (NEQ-1)

IPARAM(1) = NEQ

RPARAM(1) = HX

USER_DATA%RDATA=>RPARAM

USER_DATA%IDATA=>IPARAM

! Initial values for y, initial guesses for y'

DO I = 1, NEQ

X = (I-1) * HX

Y(I) = 1 + X

END DO

YPRIME = 0.0

NSTEPS = 10

! Always set IDO=1 on first call

IDO = 1

DO I = 1, NSTEPS

! Output solution at T=0.1,0.2,...,1.0

T = 0.1 * (I-1)

TEND = 0.1 * I

! Set IDO = 3 on last call

IF (I == NSTEPS) IDO = 3

! User-supplied Jacobian matrix (IUJAC=1)

! Banded Jacobian (MATSTR=1)

CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=1, IUJAC=1, &

MATSTR=1, ML=1, MU=1, RTOL=1.0E-4, FCN_DATA=USER_DATA)

END DO

ERRMAX = 0.0

DO I = 1, NEQ

X = (I-1) * HX

TRUE = (1+X) * EXP(T)

ERRMAX = MAX(ERRMAX, ABS(Y(I) - TRUE))

END DO

CALL UMACH(2, NOUT)

WRITE (NOUT, *) ' Max Error at T=1 is ', ERRMAX

END

SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)

USE MP_TYPES

IMPLICIT NONE

REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), HX

INTEGER IRES, LDD, I, J, NEQ, MU

TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA

NEQ = FCN_DATA%IDATA(1)

HX = FCN_DATA%RDATA(1)

MU = 1

SELECT CASE (IRES)

! F_I defined here

CASE(1)

DELTA(1) = Y(1) - EXP(T)

DO I = 2, NEQ-1

DELTA(I) = -YPRIME(I) + (Y(I+1) - 2.0 * Y(I) + Y(I-1)) &

/ HX**2 + Y(I)

END DO

DELTA(NEQ) = Y(NEQ) - 2.0 * EXP(T)

! D(I-J+MU+1,J) = D(F_I)/D(Y_J)

! in band storage mode

CASE(2)

D(MU+1,1) = 1.0

DO I = 2, NEQ-1

J = I-1

D(I-J+MU+1, J) = 1.0 / HX**2

J = I

D(I-J+MU+1, J) = -2.0 / HX**2 + 1.0

J = I+1

D(I-J+MU+1, J) = 1.0 / HX**2

END DO

D(MU+1, NEQ) = 1.0

! D(I-J+MU+1,J) = D(F_I)/D(YPRIME_J)

CASE(3)

DO I = 2, NEQ-1

D(MU+1, I) = -1.0

END DO

END SELECT

END

Output

Max Error at T=1 is 5.6743621E-5

Example 2 – Pendulum Problem

The first-order equations of motion of a point-mass m suspended on a massless wire of length L under the influence of gravity, mg, and wire tension, λ , in Cartesian coordinates (p,q) are

The problem above has an index number equal to 3, thus it cannot be solved with DAESL directly. Unfortunately, the fact that the index is greater than 1 is not obvious, but an attempt to solve it will generally produce an error message stating the corrector equation did not converge, or if IYPR=2 an error message stating that the index appears to be greater than 1 should be issued. The user then differentiates the last equation, which after replacing pʹ by u and qʹ by v, gives pu + qv = 0. This system still has index = 2 (again not obvious, the user discovers this by unsuccessfully trying to solve the new system) and the last equation must be differentiated again, to finally (after appropriate substitutions) give the equation of total energy balance:

With initial conditions and appropriate definitions of the dependent variables, the system becomes:

The initial conditions correspond to the pendulum starting in a horizontal position.

Since we have replaced the original constraint, , which requires that the pendulum length be L, by differentiating it twice, this constraint is no longer explicitly enforced, and if we try to solve the above system alone (ie, with NCON = 0), the pendulum length drifts substantially from L at larger times. DAESL therefore allows the user to add additional constraints, to be re-enforced after each time step, so we add this original constraint, as well as the intermediate constraint . Using these two supplementary constraints, (NCON = 2), the pendulum length is constant.

(Example daesl_ex2.f90)

USE DAESL_INT

USE MP_TYPES

IMPLICIT NONE

! NEQ = Number of equations

! NCON = Number of extra constraints

INTEGER, PARAMETER :: NEQ=5, NCON = 2

REAL, PARAMETER :: MASS=1.0, LENGTH=1.1, GRAVITY=9.806650

REAL T, Y(NEQ), YPRIME(NEQ), TEND, ATOL(NEQ), TOL, LEN

INTEGER NOUT, IDO, I, NSTEPS

REAL, TARGET :: RPARAM(3)

TYPE (S_FCN_DATA) USER_DATA

EXTERNAL GCN

! Pass Mass, Pendulum length and G as parameters

RPARAM(1) = MASS

RPARAM(2) = LENGTH

RPARAM(3) = GRAVITY

USER_DATA%RDATA=>RPARAM

! Initial values for y, guesses for initial y'

Y = 0.0

Y(1) = LENGTH

YPRIME = 0.0

TOL = 1.0E-5

ATOL = TOL

CALL UMACH(2, NOUT)

WRITE (NOUT, 5)

NSTEPS = 5

! Always set IDO=1 on first call

IDO = 1

DO I = 1, NSTEPS

! Output solution at T=10,20,30,40,50

T = 10.0 * (I-1)

TEND = 10.0 * I

! Set IDO = 3 on last call

IF (I.EQ.NSTEPS) IDO = 3

! User-supplied Jacobian matrix (IUJAC=1)

! Use new algorithm to get compatible y'

CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, NCON=NCON, RTOL=TOL, &

ATOL=ATOL, IYPR=2, IUJAC=1, MAXSTEPS=50000, &

FCN_DATA=USER_DATA)

! LEN = pendulum length (should be constant)

LEN = SQRT(Y(1)**2 + Y(2)**2)

WRITE (NOUT, 10) T, Y(1), Y(2), LEN

END DO

5 FORMAT (8X,'T',14X,'Y(1)',11X,'Y(2)',11X,'Length',/)

10 FORMAT (4F15.7)

END

SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)

USE MP_TYPES

IMPLICIT NONE

! Simple swinging pendulum problem

REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), MASS, &

LENGTH, GRAVITY, MG, LSQ

INTEGER IRES, LDD

TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA

MASS = FCN_DATA%RDATA(1)

LENGTH = FCN_DATA%RDATA(2)

GRAVITY = FCN_DATA%RDATA(3)

MG = MASS * GRAVITY

LSQ = LENGTH**2

SELECT CASE (IRES)

! F_I defined here

CASE(1)

DELTA(1) = Y(3) - YPRIME(1)

DELTA(2) = Y(4) - YPRIME(2)

DELTA(3) = -Y(1) * Y(5) - MASS * YPRIME(3)

DELTA(4) = -Y(2) * Y(5) - MASS * YPRIME(4) - MG

DELTA(5) = MASS * (Y(3)**2 + Y(4)**2) - MG * Y(2) - LSQ * Y(5)

! D(I,J) = D(F_I)/D(Y_J)

CASE(2)

D(1, 3) = 1.0

D(2, 4) = 1.0

D(3, 1) = -Y(5)

D(3, 5) = -Y(1)

D(4, 2) = -Y(5)

D(4, 5) = -Y(2)

D(5, 2) = -MG

D(5, 3) = MASS * 2.0 * Y(3)

D(5, 4) = MASS * 2.0 * Y(4)

D(5, 5) = -LSQ

! D(I,J) = D(F_I)/D(YPRIME_J)

CASE(3)

D(1, 1)= -1.0

D(2, 2)= -1.0

D(3, 3)= -MASS

D(4, 4)= -MASS

! DELTA(I) = D(F_I)/DT

CASE(4)

DELTA(1:5) = 0.0

! DELTA(I) = G_I

! D(I,J) = D(G_I)/D(Y_J)

CASE(5)

DELTA(1) = Y(1)**2 + Y(2)**2 - LSQ

DELTA(2) = Y(1) * Y(3) + Y(2) * Y(4)

D(1, 1) = 2.0 * Y(1)

D(1, 2) = 2.0 * Y(2)

D(1, 3) = 0.0

D(1, 4) = 0.0

D(1, 5) = 0.0

D(2, 1) = Y(3)

D(2, 2) = Y(4)

D(2, 3) = Y(1)

D(2, 4) = Y(2)

D(2, 5) = 0.0

END SELECT

END

Output

T Y(1) Y(2) Length

10.0000000 1.0998126 -0.0203017 1.0999999

20.0000000 1.0970103 -0.0810476 1.1000000

30.0000000 1.0850314 -0.1808525 1.1000004

40.0000000 1.0535675 -0.3162208 1.1000000

50.0000000 0.9896186 -0.4802662 1.1000003

Example 3 – User Solves Linear System

Consider the system of ordinary differential equations, yʹ = By, where B is the bi-diagonal matrix with (‑1, ‑1/2, ‑1/3, …, ‑1/(n‑1), 0) on the main diagonal and with 1’s along the first sub-diagonal. The initial condition is y(0) = (1,0,0,…,0)T, and since yʹ (0) = By(0) = (‑1,1,0,…,0)T, yʹ (0) is also known for this problem.

Since BT v = 0, where vi = 1/(i-1)!, v is an eigenvector of BT corresponding to the eigenvalue 0. Thus

so vT y(t) is constant. Since it has the value v T y(0) = v1 = 1 at t = 0, the constraint vTy(t) = 1 is satisfied for all t. This constraint is imposed in this example.

This example also illustrates how the user can solve his/her own linear systems (MATSTR = 2). Normally, when IRES = 6, the matrix

is computed, saved and possibly factored, using a sparse matrix factorization routine of the user’s choice. Then when IRES = 7, the system Ax = DELTA is solved, using the matrix B saved and factored earlier, and the solution is returned in DELTA. In this case, B is just a bidiagonal matrix, so there is no need to save or factor A when IRES = 6, since a bi-diagonal system can be solved directly using forward substitution, when IRES = 7.

(Example daesl_ex3.f90)

USE DAESL_INT

USE MP_TYPES

IMPLICIT NONE

! NEQ = Number of equations

INTEGER, PARAMETER :: NEQ=100

REAL T, Y(NEQ), YPRIME(NEQ), TEND, ATOL(NEQ), CON

INTEGER NOUT, IDO, I, NSTEPS

REAL, TARGET :: RPARAM(NEQ)

INTEGER, TARGET :: IPARAM(1)

TYPE (S_FCN_DATA) USER_DATA

EXTERNAL GCN

! Pass NEQ and A^T eigenvector V to GCN

IPARAM(1) = NEQ

RPARAM(1) = 1.0

DO I = 2, NEQ

RPARAM(I) = RPARAM(I-1) / FLOAT(I-1)

END DO

USER_DATA%RDATA=>RPARAM

USER_DATA%IDATA=>IPARAM

! Initial values for y, y'

Y = 0.0

Y(1) = 1.0

YPRIME = 0.0

YPRIME(1) = -1.0

YPRIME(2) = 1.0

ATOL = 1.0E-4

NSTEPS = 10

! Always set IDO=1 on first call

IDO = 1

DO I = 1, NSTEPS

! Output solution at T=1,2,...,10

T = I-1

TEND = I

! Set IDO = 3 on last call

IF (I == NSTEPS) IDO = 3

! User-defined Jacobian matrix structure (MATSTR=2)

CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=0, MATSTR=2, &

NCON=1, ATOL=ATOL, FCN_DATA=USER_DATA)

END DO

! Check if solution satisfies constraint

CON = 0.0

DO I = 1, NEQ

CON = CON + RPARAM(I) * Y(I)

END DO

CALL UMACH(2, NOUT)

WRITE (NOUT, *) ' V dot Y =', CON

END

SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)

USE MP_TYPES

IMPLICIT NONE

REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), CON, CJ

INTEGER IRES, LDD, I, NEQ

SAVE CJ

TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA

NEQ = FCN_DATA%IDATA(1)

SELECT CASE (IRES)

! F_I defined here

CASE(1)

DELTA(1) = YPRIME(1) + Y(1)

DO I = 2, NEQ-1

DELTA(I) = YPRIME(I) - Y(I-1) + Y(I) / FLOAT(I)

END DO

DELTA(NEQ) = YPRIME(NEQ) - Y(NEQ-1)

! Constraint is V dot Y = 1

CASE(5)

CON = -1.0

DO I = 1, NEQ

CON = CON + FCN_DATA%RDATA(I) * Y(I)

D(1,I) = FCN_DATA%RDATA(I)

END DO

DELTA(1) = CON

! Normally, compute matrix A = dF/dY + CJ*dF/dY'

! = -B + CJ*I here. Only CJ needs to be saved

! in this case, however, since B is bidiagonal,

! so A*x=DELTA can be solved (IRES=7) without

! saving or factoring B.

CASE(6)

CJ = DELTA(1)

! If CJ > 0 not close to zero, A is nonsingular,

! so set IRES = 0.

IF (CJ >= 1.0E-4) IRES = 0

! Solve A*x=DELTA and return x in DELTA.

CASE(7)

DELTA(1) = DELTA(1) / (1.0 + CJ)

DO I = 2, NEQ-1

DELTA(I) = (DELTA(I) + DELTA(I-1)) / (1.0 / FLOAT(I) + CJ)

END DO

DELTA(NEQ) = (DELTA(NEQ) + DELTA(NEQ-1)) / CJ

END SELECT

END

Output

V dot Y = 1.0