\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)
atomic_four_lin_ode_jac_sparsity.hpp¶
View page sourceAtomic Linear ODE Jacobian Sparsity Pattern: Example Implementation¶
Purpose¶
The jac_sparsity
routine overrides the virtual functions
used by the atomic_four base class for Jacobian sparsity calculations; see
jac_sparsity .
Notation¶
We use the notation: call_id r pattern transpose nnz , row , col , x , n , A(x) , b(x) , y(x) , m , vk(x) , and the following additional notation:
S[ g(x) ]¶
We use \(S [ g(x) ]\) to denote the sparsity pattern for a function \(g : \B{R}^n \rightarrow \B{R}^m\) as a vector of sets. To be specific, for \(i = 0, \ldots , m-1\), \(S_i [ g(x) ]\) is the set of indices between zero and \(n - 1\) such that \(\partial g_i (x) / \partial x_j\) is possibly non-zero.
N [ g(x) ]¶
We use \(N[ g(x) ]\) to denote the set of \(i\) such that \(g_i (x)\) is not identically zero.
J_i [ A(x) ]¶
We use \(J_i [ A(x) ]\) to denote the set of \(j\) between zero and \(m-1\) such that \(A_{i,j}\) is not known to be identically zero.
P_i [ g(x) ]¶
We use \(P_i [ g(x) ]\) to denote the set if sparsity pattern indices
Theory¶
This routine must calculate the following value for \(i = 0, \ldots, m-1\); see m :
The set \(S_i [ v^0 (x) ]\) has just one element corresponding to \(b_i (x)\); i.e,
see b(x) . Furthermore, for \(k > 0\),
Suppose that \(\ell\) is such that for all \(i\) the following two conditions hold
From the first condition above it follows that
Using the second condition we have
It follows that
Example¶
The file atomic_four_lin_ode_sparsity.cpp contains an example and test using this operator.
Source¶
# include <cppad/example/atomic_four/lin_ode/lin_ode.hpp>
namespace CppAD { // BEGIN_CPPAD_NAMESPACE
//
// jac_sparsity override
template <class Base>
bool atomic_lin_ode<Base>::jac_sparsity(
size_t call_id ,
bool dependency ,
const CppAD::vector<bool>& ident_zero_x ,
const CppAD::vector<bool>& select_x ,
const CppAD::vector<bool>& select_y ,
CppAD::sparse_rc< CppAD::vector<size_t> >& pattern_out )
{
//
// pattern_A, transpose, nnz
Base r;
Base step;
sparse_rc pattern_A;
bool transpose;
get(call_id, r, step, pattern_A, transpose);
size_t nnz = pattern_A.nnz();
//
// m, n
size_t m = select_y.size();
size_t n = select_x.size();
//
CPPAD_ASSERT_UNKNOWN( n == nnz + m );
CPPAD_ASSERT_UNKNOWN( pattern_A.nr() == m );
CPPAD_ASSERT_UNKNOWN( pattern_A.nc() == m );
//
// pattern_out
// Accumulates elements of the sparsity pattern for y(x) that satisfy
// select_x and select_y
pattern_out.resize(m, n, 0);
//
// list_setvec
// This vector of sets interface is not in the CppAD user API
typedef CppAD::local::sparse::list_setvec list_setvec;
//
// setvec
// Accumulates the sparsity pattern for y(x) that satisfy select_x.
// There are m sets and the set elements are between zero and n-1.
list_setvec setvec;
size_t n_set = m;
size_t end = n;
setvec.resize(n_set, end);
//
// setvec, pattern_out
// iniialize as equal to S[ v^0 (x) ]
for(size_t i = 0; i < m; ++i)
{ size_t element = nnz + i;
if( select_x[element] )
{ setvec.add_element(i, element);
if( select_y[i] )
pattern_out.push_back(i, element);
}
}
//
// non_zero
// Accumulates union_k for k < ell, N[ v^k (x) ]
// initialize as N[ v^0 (x) ]
CppAD::vector<bool> non_zero(m);
for(size_t i = 0; i < m; ++i)
non_zero[i] = ! ident_zero_x[nnz + i];
//
// change
// Did setvec or non_zero change in previous iteration of while loop
bool change = true;
while(change)
{ change = false;
// we use k = 1, 2, ... to denote the pass through this loop
// setvec[i] contains union q < k S_i [ v^q (x) ]
// non_zero contains union q < k N [ v^q (x) ]
//
// For each element of the sparsity pattern for A subject to select_x
for(size_t p = 0; p < nnz; ++p) if( select_x[p] )
{ CPPAD_ASSERT_UNKNOWN( ! ident_zero_x[p] );
size_t i = pattern_A.row()[p];
size_t j = pattern_A.col()[p];
if( transpose )
std::swap(i, j);
//
if( non_zero[j] )
{ // p, corresponding to A_{i,j}, is in P_i [ v^k (x) ]
if( ! setvec.is_element(i, p) )
{ change = true;
setvec.add_element(i, p);
if( select_y[i] )
pattern_out.push_back(i, p);
}
}
// j is in J_i [ A(x) ]
list_setvec::const_iterator itr(setvec, j);
size_t element = *itr;
while(element != end )
{ if( ! setvec.is_element(i, element) )
{ change = true;
setvec.add_element(i, element);
if( select_y[i] )
pattern_out.push_back(i, element);
}
++itr;
element = *itr;
}
// A_{i,j} update to non_zero
if( non_zero[i] == false )
{ if( non_zero[j] == true )
{ change = true;
non_zero[i] = true;
}
}
}
}
//
return true;
}
} // END_CPPAD_NAMESPACE