|
|
Prev | Next | sparse_hes_xam.m | Headings |
function ok = sparse_hes_xam()
%
% load the Cppad Swig library
m_cppad
%
% initialize return variable
ok = true;
% -----------------------------------------------------------------------
% number of dependent and independent variables
n = 3;
%
% create the independent variables ax
x = m_cppad.vec_double(n);
for i = [ 0 :(n -1) ]
x(i) = i + 2.0;
end
ax = m_cppad.independent(x);
%
% ay[i] = j * ax[j] * ax[i]
% where i = mod(j + 1, n)
ay = m_cppad.vec_a_double(n);
for j = [ 0 :(n -1) ]
i = j+1;
if( i >= n )
i = i - n;
end
aj = m_cppad.a_double(j);
ax_j = ax(j);
ax_i = ax(i);
ay(i) = aj * ax_j * ax_i;
end
%
% define af corresponding to f(x)
af = m_cppad.a_fun(ax, ay);
%
% Set select_d (domain) to all true,
% initial select_r (range) to all false
% initialize r to all zeros
select_d = m_cppad.vec_bool(n);
select_r = m_cppad.vec_bool(n);
r = m_cppad.vec_double(n);
for i = [ 0 :(n-1) ]
select_d(i) = true;
select_r(i) = false;
r(i) = 0.0;
end
%
% only select component n-1 of the range function
% f_0 (x) = (n-1) * x_{n-1} * x_0
select_r(0) = true;
r(0) = 1.0;
%
% sparisty pattern for Hessian
pattern = m_cppad.sparse_rc();
af.for_hes_sparsity(select_d, select_r, pattern);
%
% compute all possibly non-zero entries in Hessian
% (should only compute lower triangle becuase matrix is symmetric)
subset = m_cppad.sparse_rcv(pattern);
%
% work space used to save time for multiple calls
work = m_cppad.sparse_hes_work();
%
af.sparse_hes(subset, x, r, pattern, work);
%
% check that result is sparsity pattern for Hessian of f_0 (x)
ok = ok && subset.nnz() == 2 ;
row = subset.row();
col = subset.col();
val = subset.val();
for k = [ 0 :(2-1) ]
i = row(k);
j = col(k);
v = val(k);
if( i <= j )
ok = ok && i == 0;
ok = ok && j == n-1;
end
if( i >= j )
ok = ok && i == n-1;
ok = ok && j == 0;
end
ok = ok && v == n-1;
end
%
return;
end