hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_ivp_stiff_sparjac_enq (d02nr)

Purpose

nag_ode_ivp_stiff_sparjac_enq (d02nr) is an enquiry function for communicating with nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn) when supplying columns of a sparse Jacobian matrix.

Syntax

[j, iplace] = d02nr(inform)
[j, iplace] = nag_ode_ivp_stiff_sparjac_enq(inform)

Description

nag_ode_ivp_stiff_sparjac_enq (d02nr) is required when nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn) is being used with sparse matrix linear algebra. After an exit from nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn) with irevcm = 8irevcm=8, nag_ode_ivp_stiff_sparjac_enq (d02nr) must be called to determine which column of the Jacobian is required and where it is to be placed in the array rwork (a parameter of nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn)).

References

See the D02M–N sub-chapter Introduction.

Parameters

Compulsory Input Parameters

1:     inform(2323) – int64int32nag_int array
Contains information supplied by the integrator.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     j – int64int32nag_int scalar
The index jj of the column of the Jacobian which is required.
2:     iplace – int64int32nag_int scalar
Indicates which locations in the array rwork to fill with the jjth column.
If iplace = 1iplace=1, the (i,j)(i,j)th element of the Jacobian must be placed in rwork(50 + 2 × ldysav + i)rwork50+2×ldysav+i, otherwise the (i,j)(i,j)th element must be placed in rwork(50 + ldysav + i)rwork50+ldysav+i.
If jceval = 'F'jceval='F', in the previous call to nag_ode_ivp_stiff_sparjac_setup (d02nu), then iplace = 2iplace=2 always, hence the jjth column of the Jacobian must be placed in rwork(50 + ldysav + i)rwork50+ldysav+i, for i = 1,2,,neqi=1,2,,neq.

Error Indicators and Warnings

None.

Accuracy

Not applicable.

Further Comments

None.

Example

function nag_ode_ivp_stiff_sparjac_enq_example
t = 0;
tout = 10;
y = [1; 0; 0];
ydot = [0; 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');
imon = int64(0);
inln = int64(0);
ires = int64(1);
irevcm = int64(0);
lderiv = [false; false];
itask = int64(3);
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), 'Analytical', int64(100), int64(0), ...
    int64(0), int64(150), 0, 0.1, 1e-4, true, rwork);

nfails = 0;

% indexes into rwork
neq    = 3;
neqmax = neq;
lacorb = neqmax + 50;
lsavrb = lacorb + neqmax;

fprintf('\n    x          y(1)           y(2)           y(3)\n');
fprintf(' %8.3f%13.5f  %13.5f  %13.5f\n', t, y);
first_time = true;
while (irevcm ~= 0 || first_time)
  first_time = false;
  [t, tout, y, ydot, rwork, inform, ysave, wkjac, ...
   jacpvt, imon, inln, ires, irevcm, lderiv, ifail] = ...
    nag_ode_ivp_stiff_imp_revcom(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
    ysave, wkjac, jacpvt, imon, inln, ires, irevcm, lderiv, itask, itrace);

  if irevcm == 1 || irevcm == 3 || irevcm == 6 || irevcm == 11
    % Equivalent to resid evaluation in forward communication routines
    rwork(lsavrb+1) = -ydot(1) - ydot(2) - ydot(3);
    rwork(lsavrb+2) = -ydot(2);
    rwork(lsavrb+3) = -ydot(3);
    if (ires == 1)
      rwork(lsavrb+1) = rwork(lsavrb+1);
      rwork(lsavrb+2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(lsavrb+2);
      rwork(lsavrb+3) = 3e7*y(2)*y(2) + rwork(lsavrb+3);
    end
  elseif (irevcm == 2)
    % Equivalent to resid evaluation in forward communication routines
    rwork(lsavrb+1) = -rwork(lacorb+1) - rwork(lacorb+2) - rwork(lacorb+3);
    rwork(lsavrb+2) = -rwork(lacorb+2);
    rwork(lsavrb+3) = -rwork(lacorb+3);
  elseif (irevcm == 4 || irevcm == 7)
    % Equivalent to resid evaluation in forward communication routines
    rwork(lacorb+1) = -ydot(1) - ydot(2) - ydot(3);
    rwork(lacorb+2) = -ydot(2);
    rwork(lacorb+3) = -ydot(3);
    if (ires == 1)
      rwork(lacorb+1) = rwork(lacorb+1);
      rwork(lacorb+2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(lacorb+2);
      rwork(lacorb+3) = 3e7*y(2)*y(2) + rwork(lacorb+3);
    end
  elseif (irevcm == 5)
    % Equivalent to resid evaluation in forward communication routines
    ydot(1) = 0 - rwork(lsavrb+1) - rwork(lsavrb+2) - rwork(lsavrb+3);
    ydot(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) - rwork(lsavrb+2);
    ydot(3) = 3e7*y(2)*y(2) - rwork(lsavrb+3);
  elseif (irevcm == 8)
    % Equivalent to jac evaluation in forward communication routines
    [j, iplace] = nag_ode_ivp_stiff_sparjac_enq(inform);
    hxd = rwork(16)*rwork(20);
    if (iplace < 2)
      if (j < 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 0 - hxd*(0.04);
        rwork(lsavrb+3) = 0 - hxd*(0);
      elseif (j == 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 1 - hxd*(-1e4*y(3)-6e7*y(2));
        rwork(lsavrb+3) = 0 - hxd*(6e7*y(2));
      elseif (j > 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 0 - hxd*(-1e4*y(2));
        rwork(lsavrb+3) = 1 - hxd*(0);
      end
    else
      if (j < 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 0 - hxd*(0.04);
        rwork(lacorb+3) = 0 - hxd*(0);
      elseif (j == 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 1 - hxd*(-1e4*y(3)-6e7*y(2));
        rwork(lacorb+3) = 0 - hxd*(6e7*y(2));
      elseif (j > 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 0 - hxd*(-1e4*y(2));
        rwork(lacorb+3) = 1 - hxd*(0);
      end
    end
  elseif (irevcm == 10)
    % Step failure
    nfails = nfails + 1;
  end
end

if (ifail == 0)
  [hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
    nag_ode_ivp_stiff_integ_diag(int64(neq), int64(neqmax), rwork, inform);

  [y, ifail] = nag_ode_ivp_stiff_interp(tout, int64(neq), int64(neq), ysave, rwork);

  fprintf(' %8.3f%13.5f  %13.5f  %13.5f\n', tout, y);
  fprintf('\nhused = %12.5e hnext = %12.5e tcur = %12.5e\n', hu, h, tcur);
  fprintf('nst = %6d   nre = %6d   nje = %6d\n', nst, nre, nje);
  fprintf('nqa = %6d   nq  = %6d niter = %6d\n', nqu, nq, niter);
  fprintf('Max err comp = %4d   No. of failed steps = %4d\n', imxer, nfails);

  icall = int64(0);
  [liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock] = ...
    nag_ode_ivp_stiff_sparjac_diag(icall, true, inform);

  fprintf('\nnjcpvt (required %4d  used %8d)\n', liwreq, liwusd);
  fprintf('nwjac (required %4d  used %8d)\n', lrwreq, lrwusd);
  fprintf('No. of LU-decomps %4d  No. of nonzeros %8d\n', nlu, nz)
  fprintf('No. of fcn calls to form Jacobian %4d Try isplit %4d\n', ngp, isplit);
  fprintf('Growth est %8d  No. of blocks on diagonal %4d\n', igrow, nblock);
elseif (ifail == 10)
  icall = int64(1)
  [liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock] = ...
    nag_ode_ivp_stiff_sparjac_diag(icall, true, inform);

  fprintf('\nnjcpvt (required %4d  used %8d)\n', liwreq, liwusd);
  fprintf('nwjac (required %4d  used %8d)\n', lrwreq, lrwusd);
else
  fprintf('Exit nag_ode_ivp_stiff_imp_revcom with ifail = %d and t = %12.5e\n', ifail, t);
end
 

    x          y(1)           y(2)           y(3)
    0.000      1.00000        0.00000        0.00000
 WARNING... EQUATION(=I1) AND POSSIBLY OTHER EQUATIONS ARE
 IMPLICIT AND IN CALCULATING THE INITIAL VALUES THE EQNS
 WILL BE TREATED AS IMPLICIT.
 IN ABOVE MESSAGE I1 =         1
   10.000      0.84135        0.00002        0.15863

hused =  9.01811e-01 hnext =  9.01811e-01 tcur =  1.07667e+01
nst =     55   nre =     89   nje =     17
nqa =      4   nq  =      4 niter =     80
Max err comp =    3   No. of failed steps =    4

njcpvt (required   99  used      150)
nwjac (required   31  used       75)
No. of LU-decomps   17  No. of nonzeros        8
No. of fcn calls to form Jacobian    0 Try isplit   73
Growth est     1531  No. of blocks on diagonal    1

function d02nr_example
t = 0;
tout = 10;
y = [1; 0; 0];
ydot = [0; 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');
imon = int64(0);
inln = int64(0);
ires = int64(1);
irevcm = int64(0);
lderiv = [false; false];
itask = int64(3);
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), 'Analytical', int64(100), int64(0), ...
    int64(0), int64(150), 0, 0.1, 1e-4, true, rwork);

nfails = 0;

% indexes into rwork
neq    = 3;
neqmax = neq;
lacorb = neqmax + 50;
lsavrb = lacorb + neqmax;

fprintf('\n    x          y(1)           y(2)           y(3)\n');
fprintf(' %8.3f%13.5f  %13.5f  %13.5f\n', t, y);
first_time = true;
while (irevcm ~= 0 || first_time)
  first_time = false;
  [t, tout, y, ydot, rwork, inform, ysave, wkjac, ...
   jacpvt, imon, inln, ires, irevcm, lderiv, ifail] = ...
    d02nn(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
    ysave, wkjac, jacpvt, imon, inln, ires, irevcm, lderiv, itask, itrace);

  if irevcm == 1 || irevcm == 3 || irevcm == 6 || irevcm == 11
    % Equivalent to resid evaluation in forward communication routines
    rwork(lsavrb+1) = -ydot(1) - ydot(2) - ydot(3);
    rwork(lsavrb+2) = -ydot(2);
    rwork(lsavrb+3) = -ydot(3);
    if (ires == 1)
      rwork(lsavrb+1) = rwork(lsavrb+1);
      rwork(lsavrb+2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(lsavrb+2);
      rwork(lsavrb+3) = 3e7*y(2)*y(2) + rwork(lsavrb+3);
    end
  elseif (irevcm == 2)
    % Equivalent to resid evaluation in forward communication routines
    rwork(lsavrb+1) = -rwork(lacorb+1) - rwork(lacorb+2) - rwork(lacorb+3);
    rwork(lsavrb+2) = -rwork(lacorb+2);
    rwork(lsavrb+3) = -rwork(lacorb+3);
  elseif (irevcm == 4 || irevcm == 7)
    % Equivalent to resid evaluation in forward communication routines
    rwork(lacorb+1) = -ydot(1) - ydot(2) - ydot(3);
    rwork(lacorb+2) = -ydot(2);
    rwork(lacorb+3) = -ydot(3);
    if (ires == 1)
      rwork(lacorb+1) = rwork(lacorb+1);
      rwork(lacorb+2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) + rwork(lacorb+2);
      rwork(lacorb+3) = 3e7*y(2)*y(2) + rwork(lacorb+3);
    end
  elseif (irevcm == 5)
    % Equivalent to resid evaluation in forward communication routines
    ydot(1) = 0 - rwork(lsavrb+1) - rwork(lsavrb+2) - rwork(lsavrb+3);
    ydot(2) = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)*y(2) - rwork(lsavrb+2);
    ydot(3) = 3e7*y(2)*y(2) - rwork(lsavrb+3);
  elseif (irevcm == 8)
    % Equivalent to jac evaluation in forward communication routines
    [j, iplace] = d02nr(inform);
    hxd = rwork(16)*rwork(20);
    if (iplace < 2)
      if (j < 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 0 - hxd*(0.04);
        rwork(lsavrb+3) = 0 - hxd*(0);
      elseif (j == 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 1 - hxd*(-1e4*y(3)-6e7*y(2));
        rwork(lsavrb+3) = 0 - hxd*(6e7*y(2));
      elseif (j > 2)
        rwork(lsavrb+1) = 1 - hxd*(0);
        rwork(lsavrb+2) = 0 - hxd*(-1e4*y(2));
        rwork(lsavrb+3) = 1 - hxd*(0);
      end
    else
      if (j < 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 0 - hxd*(0.04);
        rwork(lacorb+3) = 0 - hxd*(0);
      elseif (j == 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 1 - hxd*(-1e4*y(3)-6e7*y(2));
        rwork(lacorb+3) = 0 - hxd*(6e7*y(2));
      elseif (j > 2)
        rwork(lacorb+1) = 1 - hxd*(0);
        rwork(lacorb+2) = 0 - hxd*(-1e4*y(2));
        rwork(lacorb+3) = 1 - hxd*(0);
      end
    end
  elseif (irevcm == 10)
    % Step failure
    nfails = nfails + 1;
  end
end

if (ifail == 0)
  [hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
    d02ny(int64(neq), int64(neqmax), rwork, inform);

  [y, ifail] = d02mz(tout, int64(neq), int64(neq), ysave, rwork);

  fprintf(' %8.3f%13.5f  %13.5f  %13.5f\n', tout, y);
  fprintf('\nhused = %12.5e hnext = %12.5e tcur = %12.5e\n', hu, h, tcur);
  fprintf('nst = %6d   nre = %6d   nje = %6d\n', nst, nre, nje);
  fprintf('nqa = %6d   nq  = %6d niter = %6d\n', nqu, nq, niter);
  fprintf('Max err comp = %4d   No. of failed steps = %4d\n', imxer, nfails);

  icall = int64(0);
  [liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock] = ...
    d02nx(icall, true, inform);

  fprintf('\nnjcpvt (required %4d  used %8d)\n', liwreq, liwusd);
  fprintf('nwjac (required %4d  used %8d)\n', lrwreq, lrwusd);
  fprintf('No. of LU-decomps %4d  No. of nonzeros %8d\n', nlu, nz)
  fprintf('No. of fcn calls to form Jacobian %4d Try isplit %4d\n', ngp, isplit);
  fprintf('Growth est %8d  No. of blocks on diagonal %4d\n', igrow, nblock);
elseif (ifail == 10)
  icall = int64(1)
  [liwreq, liwusd, lrwreq, lrwusd, nlu, nz, ngp, isplit, igrow, nblock] = ...
    d02nx(icall, true, inform);

  fprintf('\nnjcpvt (required %4d  used %8d)\n', liwreq, liwusd);
  fprintf('nwjac (required %4d  used %8d)\n', lrwreq, lrwusd);
else
  fprintf('Exit d02nn with ifail = %d and t = %12.5e\n', ifail, t);
end
 

    x          y(1)           y(2)           y(3)
    0.000      1.00000        0.00000        0.00000
 WARNING... EQUATION(=I1) AND POSSIBLY OTHER EQUATIONS ARE
 IMPLICIT AND IN CALCULATING THE INITIAL VALUES THE EQNS
 WILL BE TREATED AS IMPLICIT.
 IN ABOVE MESSAGE I1 =         1
   10.000      0.84135        0.00002        0.15863

hused =  9.01811e-01 hnext =  9.01811e-01 tcur =  1.07667e+01
nst =     55   nre =     89   nje =     17
nqa =      4   nq  =      4 niter =     80
Max err comp =    3   No. of failed steps =    4

njcpvt (required   99  used      150)
nwjac (required   31  used       75)
No. of LU-decomps   17  No. of nonzeros        8
No. of fcn calls to form Jacobian    0 Try isplit   73
Growth est     1531  No. of blocks on diagonal    1


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013