\(\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_for_type.hpp¶
View page sourceAtomic Linear ODE Forward Type Calculation: Example Implementation¶
Purpose¶
The for_type
routine overrides the virtual functions
used by the atomic_four base; see
for_type .
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:
T(s)¶
We use \(\R{T} ( s )\) to denote the ad_type of a scalar value \(s\). There are four possible ad_types : identical_zero, constant, dynamic, and variable in that order.
Theory¶
This routine must calculate the following value for \(i = 0, \ldots, m-1\); see m :
The type \(\R{T} [ v_i^0 (x) ] = \R{T}[ b_i (x) ]\). This is easy to calculate given the type of the components of x ; see b(x) . Furthermore, for \(k > 0\)
If \(A_{i,j} (x)\) is not in the sparsity pattern , it is identically zero. Furthermore we are allowing for the case where \(A_{i,j} (x)\) is in the pattern and it is identically zero; i.e., the sparsity pattern is not efficient as it could be. The type \(\R{T} [ A_{i,j} (x) ]\) for components in the sparsity pattern is easy to calculate given the type of the components of x ; see A(x) . Suppose \(\ell\) is such that for all \(i\)
It follows that
From this it is clear that
Source¶
# include <cppad/example/atomic_four/lin_ode/lin_ode.hpp>
namespace CppAD { // BEGIN_CPPAD_NAMESPACE
//
// for_type override
template <class Base>
bool atomic_lin_ode<Base>::for_type(
size_t call_id ,
const CppAD::vector<CppAD::ad_type_enum>& type_x ,
CppAD::vector<CppAD::ad_type_enum>& type_y )
{
// pattern, transpose, nnz
Base r;
Base step;
sparse_rc pattern;
bool transpose;
get(call_id, step, r, pattern, transpose);
size_t nnz = pattern.nnz();
//
// m
size_t m = type_y.size();
CPPAD_ASSERT_UNKNOWN( pattern.nr() == m );
CPPAD_ASSERT_UNKNOWN( pattern.nc() == m );
//
// type_x
CPPAD_ASSERT_UNKNOWN( type_x.size() == nnz + m );
//
// type_y[i] = type_b[i]
// type_y[i] = max T[ v_i^k (x) ] for k = 0
for(size_t i = 0; i < m; ++i)
type_y[i] = type_x[nnz + i];
//
// change
// Did type_y change during the previous iteration of the while loop
bool change = true;
while(change)
{ change = false;
// we use k = 1, 2, ... to denote the pass through this loop
// type_y[i] = max q < k T[ v_i^q (x) ]
//
for(size_t p = 0; p < nnz; ++p)
{ size_t i = pattern.row()[p];
size_t j = pattern.col()[p];
if( transpose )
std::swap(i, j);
//
// type_y[i], change
if( type_x[p] != identical_zero_enum )
{ // A_ij (x) is not identically zero
if( type_y[j] > type_y[i] )
{ change = true;
type_y[i] = type_y[j];
}
}
//
// type_y[i], change
if( type_y[j] != identical_zero_enum )
{ // There is a q < k such that v_j^q (x) not identically zero
if( type_x[p] > type_y[i] )
{ change = true;
type_y[i] = type_x[p];
}
}
}
}
return true;
}
} // END_CPPAD_NAMESPACE