\(\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_hes_sparsity.hpp¶
View page sourceAtomic Linear ODE Hessian Sparsity Pattern: Example Implementation¶
Purpose¶
The hes_sparsity
routine overrides the virtual functions
used by the atomic_four base class for Hessian sparsity calculations; see
hes_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:
wk(x)¶
Because we are using the Rosen34 solver, our actual sequence of operations is only fourth order accurate. So it suffices to compute the sparsity pattern for
Note that the factor \(r / k\), in the definition of \(v^k (x)\), is constant (with respect to the variables). Hence it suffices to compute the sparsity pattern for
where \(w^0 (x) = b(x)\) and for \(k = 1, 2, \ldots\), \(w^k (x) = A(x) w^{k-1} (x)\).
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
//
// hes_sparsity override
template <class Base>
bool atomic_lin_ode<Base>::hes_sparsity(
size_t call_id ,
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 )
{ //
// adouble
typedef AD<double> adouble;
//
// 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();
//
// au
vector<adouble> au(n);
for(size_t j = 0; j < n; ++j)
au[j] = 1.0;
Independent( au );
//
// ax
vector<adouble> ax(n);
for(size_t j = 0; j < n; ++j)
if( ident_zero_x[j] )
ax[j] = 0.0;
else
ax[j] = au[j];
//
// aw
// aw = w^0 (x)
vector<adouble> aw(m);
for(size_t i = 0; i < m; ++i)
aw[i] = ax[nnz + i];
//
// ah = w^0 (x)
vector<adouble> ah = aw;
//
// ah = sum_k=0^4 w^k (x)
vector<adouble> awk(m);
for(size_t k = 1; k < 5; ++k)
{ // aw = w^{k-1} (x)
//
// awk = w^k (x)
for(size_t i = 0; i < m; ++i)
awk[i] = 0.0;
for(size_t p = 0; p < nnz; ++p)
{ //
// i, j
size_t i = pattern_A.row()[p];
size_t j = pattern_A.col()[p];
if( transpose )
std::swap(i, j);
//
// awk
awk[i] += ax[p] * aw[j];
}
//
// ah = ah + w^k(x)
for(size_t i = 0; i < m; ++i)
ah[i] += awk[i];
//
// aw = w^k (x)
aw = awk;
}
//
// h
ADFun<double> h;
h.Dependent(au, ah);
//
// pattern_out
// can use h.for_hes_sparsity or h.rev_hes_sparsity
bool internal_bool = false;
h.for_hes_sparsity(select_x, select_y, internal_bool, pattern_out);
//
return true;
}
} // END_CPPAD_NAMESPACE