Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_ivp_stiff_sparjac_setup (d02nu)

## Purpose

nag_ode_ivp_stiff_sparjac_setup (d02nu) is a setup function which must be called prior to an integrator in sub-chapter D02M–N, if sparse matrix linear algebra is required.

## Syntax

[jacpvt, rwork, ifail] = d02nu(neq, neqmax, jceval, nwkjac, ia, ja, njcpvt, sens, u, eta, lblock, rwork, 'nia', nia, 'nja', nja, 'isplit', isplit)
[jacpvt, rwork, ifail] = nag_ode_ivp_stiff_sparjac_setup(neq, neqmax, jceval, nwkjac, ia, ja, njcpvt, sens, u, eta, lblock, rwork, 'nia', nia, 'nja', nja, 'isplit', isplit)

## Description

nag_ode_ivp_stiff_sparjac_setup (d02nu) defines the linear algebra to be used as sparse matrix linear algebra, permits you to specify the method for calculating the Jacobian and its structure, and checks the validity of certain input values.

## References

See the D02M–N sub-chapter Introduction.

## Parameters

### Compulsory Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential equations.
Constraint: 1neqneqmax$1\le {\mathbf{neq}}\le {\mathbf{neqmax}}$.
2:     neqmax – int64int32nag_int scalar
A bound on the maximum number of differential equations to be solved during the integration.
Constraint: ${\mathbf{neqmax}}\ge {\mathbf{neq}}$.
3:     jceval – string (length ≥ 1)
Specifies the technique to be used to compute the Jacobian.
jceval = 'N'${\mathbf{jceval}}=\text{'N'}$
The sparsity structure and the value of the Jacobian are to be determined numerically by the integrator.
jceval = 'S'${\mathbf{jceval}}=\text{'S'}$
The sparsity structure of the Jacobian is supplied in the arrays ia and ja but its value is to be determined numerically. This is the recommended mode of operation unless it is a simple matter to supply the Jacobian.
jceval = 'A'${\mathbf{jceval}}=\text{'A'}$
The Jacobian will be evaluated by calls to jac. The sparsity structure will be estimated by calls to jac; that is, no explicit sparsity structure need be supplied in the arrays ia and ja.
jceval = 'F'${\mathbf{jceval}}=\text{'F'}$
The sparsity structure of the Jacobian is supplied in ia and ja, and its value will be determined by calls to jac. This is the recommended mode of operation if the jac is simple to form.
jceval = 'D'${\mathbf{jceval}}=\text{'D'}$
The default choice is to be made. In this case 'D' is interpreted as 'S'.
If the sparsity structure is supplied in arrays ia and ja, then any evidence from the numerical or analytical formation of the Jacobian that this structure is not correct, is ignored.
Only the first character of the actual parameter jceval is passed to nag_ode_ivp_stiff_sparjac_setup (d02nu); hence it is permissible for the actual argument to be more descriptive, e.g., ‘Numerical’, ‘Structural’, ‘Analytical’, ‘Full information’ or ‘Default’ in a call to nag_ode_ivp_stiff_sparjac_setup (d02nu).
If the option jceval = 'N'${\mathbf{jceval}}=\text{'N'}$, 'S'$\text{'S'}$ or 'D'$\text{'D'}$ is used then the actual argument corresponding to jac in the call to nag_ode_ivp_stiff_exp_sparjac (d02nd) or nag_ode_ivp_stiff_imp_sparjac (d02nj) must be either nag_ode_ivp_stiff_exp_sparjac_dummy_jac (d02ndz) or nag_ode_ivp_stiff_imp_sparjac_dummy_jac (d02njz) respectively.
If integration is to be performed by reverse communication (nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn)) then jceval should be set to either 'N' or 'A'. In this case ia and ja are not used and their lengths may be set to 1$1$.
Constraint: jceval = 'N'${\mathbf{jceval}}=\text{'N'}$, 'S'$\text{'S'}$, 'A'$\text{'A'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$.
4:     nwkjac – int64int32nag_int scalar
The size of the array wkjac, which you are supplying to the integrator, as declared in the (sub)program from which nag_ode_ivp_stiff_sparjac_setup (d02nu) is called.
Suggested value: nwkjac = 4 × neqmax${\mathbf{nwkjac}}=4×{\mathbf{neqmax}}$ if jceval = 'N'${\mathbf{jceval}}=\text{'N'}$ or 'A'$\text{'A'}$. If nwkjac is less than this estimate, then a message is printed on the current advisory message unit (see nag_file_set_unit_advisory (x04ab)), and execution continues.
Constraint: if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, nwkjacnelement + 2 × neq${\mathbf{nwkjac}}\ge \mathit{nelement}+2×{\mathbf{neq}}$, where nelement$\mathit{nelement}$ is the total number of nonzeros.
5:     ia(nia) – int64int32nag_int array
nia, the dimension of the array, must satisfy the constraint
• if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, nianeq + 1${\mathbf{nia}}\ge {\mathbf{neq}}+1$;
• otherwise nia1${\mathbf{nia}}\ge 1$.
If jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, ia must contain details of the sparsity pattern to be used for the Jacobian. See ja.
ia is not used if jceval = 'N'${\mathbf{jceval}}=\text{'N'}$ or 'A'$\text{'A'}$.
6:     ja(nja) – int64int32nag_int array
nja, the dimension of the array, must satisfy the constraint
• if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, njaia(neq + 1)1${\mathbf{nja}}\ge {\mathbf{ia}}\left({\mathbf{neq}}+1\right)-1$;
• otherwise nja1${\mathbf{nja}}\ge 1$.
If jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, ja must contain details of the sparsity pattern to be used for the Jacobian. ja contains the row indices where nonzero elements occur, reading in column-wise order, and ia contains the starting locations in ja of the descriptions of columns 1,2,,neq$1,2,\dots ,{\mathbf{neq}}$ in that order, with ia(1) = 1${\mathbf{ia}}\left(1\right)=1$. Thus for each column index j = 1,2,,neq$j=1,2,\dots ,{\mathbf{neq}}$, the values of the row index i$i$ in column j$j$ where a nonzero element may occur are given by
 i = ja(k) $i=jak$
where ia(j)k < ia(j + 1)${\mathbf{ia}}\left(j\right)\le k<{\mathbf{ia}}\left(j+1\right)$.
Thus the total number of nonzeros, nelement$\mathit{nelement}$, must be ia(neq + 1)1${\mathbf{ia}}\left({\mathbf{neq}}+1\right)-1$. For example, for the following matrix
 x 0 x 0 0 0 x x x 0 x x x 0 0 x 0 0 x x 0 0 0 x x
$x 0 x 0 0 0 x x x 0 x x x 0 0 x 0 0 x x 0 0 0 x x$
where x$x$ represents nonzero elements (13 in all) the arrays ia and ja should be
 ia(k) 1 4 6 9 12 14 ja(k) 1 3 4 2 03 01 2 3 2 4 5 4 5
$iak 1 4 6 9 12 14 jak 1 3 4 2 03 01 2 3 2 4 5 4 5$
ja is not used if jceval = 'N'${\mathbf{jceval}}=\text{'N'}$ or 'A'$\text{'A'}$.
7:     njcpvt – int64int32nag_int scalar
The length of the array jacpvt, which you are supplying to the integrator, as dimensioned in the sub(program) from which nag_ode_ivp_stiff_sparjac_setup (d02nu) is called.
Suggested value: njcpvt = 20 × neqmax${\mathbf{njcpvt}}=20×{\mathbf{neqmax}}$ if jceval = 'N'${\mathbf{jceval}}=\text{'N'}$ or 'A'$\text{'A'}$. If njcpvt is less than this estimate, then a message is printed on the current advisory message unit (see nag_file_set_unit_advisory (x04ab)), and execution continues.
Constraint: if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, njcpvt3 × nelement + 14 × neq${\mathbf{njcpvt}}\ge 3×\mathit{nelement}+14×{\mathbf{neq}}$, where nelement$\mathit{nelement}$ is the total number of nonzeros.
8:     sens – double scalar
A threshold parameter used to determine whether or not a matrix element is zero; when sens is set to 0.0$0.0$ on entry, the function will use sens = 100.0 × machine precision. Otherwise the absolute value of sens is used.
9:     u – double scalar
Should have a value between 0.0$0.0$ and 0.9999$0.9999$. Otherwise a default value of 0.1$0.1$ is used. When the sparsity pattern has been evaluated, the first Jacobian computed is decomposed with u governing the choice of pivots; subsequent Jacobian decompositions use the same pattern of decomposition until the sparsity pattern is re-evaluated. When searching a row for a pivot, any element is excluded from the search which is less than u times the largest of those elements in the row available as pivots. Thus decreasing u biases the algorithm towards maintaining sparsity at the expense of numerical stability.
10:   eta – double scalar
A relative pivot threshold, below which on subsequent decompositions (as described under u), an internal error is provoked.
eta > 1.0${\mathbf{eta}}>1.0$
No check on pivot size is made.
eta0.0${\mathbf{eta}}\le 0.0$
The default value eta = 1.0e−4${\mathbf{eta}}=\text{1.0e−4}$ is used.
11:   lblock – logical scalar
Indicates if preordering is used before decomposition.
If lblock = true${\mathbf{lblock}}=\mathbf{true}$, on entry, the Jacobian matrix is preordered to block lower triangular form before a decomposition is performed (this is the recommended mode). If you know the structure of the Jacobian to be irreducible, that is not permutable to block lower triangular form, then you should set lblock = false${\mathbf{lblock}}=\mathbf{false}$. For example, a Jacobian arising from using the method of lines for parabolic partial differential equations would normally be irreducible. (See the specification of nag_ode_ivp_stiff_sparjac_diag (d02nx) for optional output concerning lblock.)
12:   rwork(50 + 4 × neqmax$50+4×{\mathbf{neqmax}}$) – double array
This must be the same workspace array as the array rwork supplied to the integrator. It is used to pass information from the setup function to the integrator and therefore the contents of this array must not be changed before calling the integrator.

### Optional Input Parameters

1:     nia – int64int32nag_int scalar
Default: The dimension of the array ia.
The dimension of the array ia as declared in the (sub)program from which nag_ode_ivp_stiff_sparjac_setup (d02nu) is called.
Constraints:
• if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, nianeq + 1${\mathbf{nia}}\ge {\mathbf{neq}}+1$;
• otherwise nia1${\mathbf{nia}}\ge 1$.
2:     nja – int64int32nag_int scalar
Default: The dimension of the array ja.
The dimension of the array ja as declared in the (sub)program from which nag_ode_ivp_stiff_sparjac_setup (d02nu) is called.
Constraints:
• if jceval = 'S'${\mathbf{jceval}}=\text{'S'}$, 'F'$\text{'F'}$ or 'D'$\text{'D'}$, njaia(neq + 1)1${\mathbf{nja}}\ge {\mathbf{ia}}\left({\mathbf{neq}}+1\right)-1$;
• otherwise nja1${\mathbf{nja}}\ge 1$.
3:     isplit – int64int32nag_int scalar
This parameter is used for splitting the integer workspace jacpvt to effect an efficient decomposition. It must satisfy 1isplit99$1\le {\mathbf{isplit}}\le 99$. If isplit lies outside this range on entry, a default value of 73$73$ is used. An appropriate value for isplit for subsequent runs on similar problems is available via the optional output nag_ode_ivp_stiff_sparjac_diag (d02nx).
Default: 73$73$

None.

### Output Parameters

1:     jacpvt(njcpvt) – int64int32nag_int array
Data relating to the Jacobian sparsity structure.
2:     rwork(50 + 4 × neqmax$50+4×{\mathbf{neqmax}}$) – double array
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
On entry, an illegal input was detected.

## Accuracy

Not applicable.

nag_ode_ivp_stiff_sparjac_setup (d02nu) must be called as a setup function before a call to either nag_ode_ivp_stiff_exp_sparjac (d02nd) or nag_ode_ivp_stiff_imp_sparjac (d02nj) and may be called as the linear algebra setup function before a call to nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn).

## Example

function nag_ode_ivp_stiff_sparjac_setup_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62, 1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(100, 1);
jacpvt = zeros(150, 1, 'int64');
itrace = int64(0);
[const, rwork, ifail] = ...
nag_ode_ivp_stiff_bdf(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
0, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);
[jacpvt, rwork, ifail] = ...
nag_ode_ivp_stiff_sparjac_setup(int64(3), int64(3), 'Numerical', int64(100), int64(0), ...
int64(0), int64(150), 0, 0.1, 1e-4, true, rwork);
[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ...
jacpvtOut, ifail] = ...
nag_ode_ivp_stiff_exp_sparjac(t, tout, y, rwork, rtol, atol, itol, inform, ...
@fcn, ysave, 'nag_ode_ivp_stiff_exp_sparjac_dummy_jac', wkjac, jacpvt, 'nag_ode_ivp_stiff_exp_fulljac_dummy_monit', ...
tOut, yOut, ydot, informOut, ifail

function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);

tOut =

10

yOut =

0.8414
0.0000
0.1586

ydot =

-0.0075
-0.0000
0.0075

informOut =

55
136
16
4
4
78
0
3
0
100
150
29
71
16
7
3
73
1108
1
0
0
0
0

ifail =

0

function d02nu_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62, 1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(100, 1);
jacpvt = zeros(150, 1, 'int64');
itrace = int64(0);
[const, rwork, ifail] = ...
d02nv(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
0, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);
[jacpvt, rwork, ifail] = ...
d02nu(int64(3), int64(3), 'Numerical', int64(100), int64(0), ...
int64(0), int64(150), 0, 0.1, 1e-4, true, rwork);
[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ...
jacpvtOut, ifail] = ...
d02nd(t, tout, y, rwork, rtol, atol, itol, inform, ...
@fcn, ysave, 'd02ndz', wkjac, jacpvt, 'd02nby', ...
tOut, yOut, ydot, informOut, ifail

function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);

tOut =

10

yOut =

0.8414
0.0000
0.1586

ydot =

-0.0075
-0.0000
0.0075

informOut =

55
136
16
4
4
78
0
3
0
100
150
29
71
16
7
3
73
1108
1
0
0
0
0

ifail =

0