Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_eigen_mat_mul.cpp ADValued-> Arithmetic unary_standard_math binary_math CondExp Discrete numeric_limits atomic atomic-> checkpoint atomic_base atomic_base-> atomic_ctor atomic_option atomic_afun atomic_forward atomic_reverse atomic_for_sparse_jac atomic_rev_sparse_jac atomic_for_sparse_hes atomic_rev_sparse_hes atomic_base_clear atomic_get_started.cpp atomic_norm_sq.cpp atomic_reciprocal.cpp atomic_set_sparsity.cpp atomic_tangent.cpp atomic_eigen_mat_mul.cpp atomic_eigen_mat_inv.cpp atomic_eigen_cholesky.cpp atomic_mat_mul.cpp atomic_eigen_mat_mul.cpp-> atomic_eigen_mat_mul.hpp Headings-> Description Class Definition Use Atomic Function ---..Constructor

$\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 Eigen Matrix Multiply: Example and Test

Description
The ADFun function object f for this example is $$f(x) = \left( \begin{array}{cc} 0 & 0 \\ 1 & 2 \\ x_0 & x_1 \end{array} \right) \left( \begin{array}{c} x_0 \\ x_1 \end{array} \right) = \left( \begin{array}{c} 0 \\ x_0 + 2 x_1 \\ x_0 x_0 + x_1 x_1 ) \end{array} \right)$$

Class Definition
This example uses the file atomic_eigen_mat_mul.hpp which defines matrix multiply as a atomic_base operation.

Use Atomic Function
# include <cppad/cppad.hpp>

bool eigen_mat_mul(void)
{     //
typedef double                                            scalar;
//
bool ok    = true;
scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
//

Constructor
     // -------------------------------------------------------------------
// object that multiplies arbitrary matrices
atomic_eigen_mat_mul<scalar> mat_mul;
// -------------------------------------------------------------------
// declare independent variable vector x
size_t n = 2;
for(size_t j = 0; j < n; j++)
// -------------------------------------------------------------------
//        [ 0     0    ]
// left = [ 1     2    ]
//        [ x[0]  x[1] ]
size_t nr_left  = 3;
size_t n_middle   = 2;
// -------------------------------------------------------------------
// right = [ x[0] , x[1] ]^T
size_t nc_right = 1;
// -------------------------------------------------------------------
// use atomic operation to multiply left * right
// -------------------------------------------------------------------
// check that first component of result is a parameter
// and the other components are varaibles.
ok &= Parameter( ad_result(0, 0) );
ok &= Variable(  ad_result(1, 0) );
ok &= Variable(  ad_result(2, 0) );
// -------------------------------------------------------------------
// declare the dependent variable vector y
size_t m = 3;
for(size_t i = 0; i < m; i++)
// -------------------------------------------------------------------
// check zero order forward mode
for(size_t i = 0; i < n; i++)
x[i] = scalar(i + 2);
y   = f.Forward(0, x);
ok &= NearEqual(y[0], 0.0,                       eps, eps);
ok &= NearEqual(y[1], x[0] + 2.0 * x[1],         eps, eps);
ok &= NearEqual(y[2], x[0] * x[0] + x[1] * x[1], eps, eps);
// -------------------------------------------------------------------
// check first order forward mode
x1[0] = 1.0;
x1[1] = 0.0;
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0,        eps, eps);
ok   &= NearEqual(y1[1], 1.0,        eps, eps);
ok   &= NearEqual(y1[2], 2.0 * x[0], eps, eps);
x1[0] = 0.0;
x1[1] = 1.0;
y1    = f.Forward(1, x1);
ok   &= NearEqual(y1[0], 0.0,        eps, eps);
ok   &= NearEqual(y1[1], 2.0,        eps, eps);
ok   &= NearEqual(y1[2], 2.0 * x[1], eps, eps);
// -------------------------------------------------------------------
// check second order forward mode
x2[0] = 0.0;
x2[1] = 0.0;
y2    = f.Forward(2, x2);
ok   &= NearEqual(y2[0], 0.0, eps, eps);
ok   &= NearEqual(y2[1], 0.0, eps, eps);
ok   &= NearEqual(y2[2], 1.0, eps, eps); // 1/2 * f_1''(x)
// -------------------------------------------------------------------
// check first order reverse mode
w[0]  = 0.0;
w[1]  = 1.0;
w[2]  = 0.0;
d1w   = f.Reverse(1, w);
ok   &= NearEqual(d1w[0], 1.0, eps, eps);
ok   &= NearEqual(d1w[1], 2.0, eps, eps);
w[0]  = 0.0;
w[1]  = 0.0;
w[2]  = 1.0;
d1w   = f.Reverse(1, w);
ok   &= NearEqual(d1w[0], 2.0 * x[0], eps, eps);
ok   &= NearEqual(d1w[1], 2.0 * x[1], eps, eps);
// -------------------------------------------------------------------
// check second order reverse mode
d2w   = f.Reverse(2, w);
// partial f_2 w.r.t. x_0
ok   &= NearEqual(d2w[0 * 2 + 0], 2.0 * x[0], eps, eps);
// partial f_2 w.r.t  x_1
ok   &= NearEqual(d2w[1 * 2 + 0], 2.0 * x[1], eps, eps);
// partial f_2 w.r.t x_1, x_0
ok   &= NearEqual(d2w[0 * 2 + 1], 0.0,        eps, eps);
// partial f_2 w.r.t x_1, x_1
ok   &= NearEqual(d2w[1 * 2 + 1], 2.0,        eps, eps);
// -------------------------------------------------------------------
// check forward Jacobian sparsity
std::set<size_t> check_set;
for(size_t j = 0; j < n; j++)
r[j].insert(j);
s      = f.ForSparseJac(n, r);
check_set.clear();
ok    &= s[0] == check_set;
check_set.insert(0);
check_set.insert(1);
ok    &= s[1] == check_set;
ok    &= s[2] == check_set;
// -------------------------------------------------------------------
// check reverse Jacobian sparsity
r.resize(m);
for(size_t i = 0; i < m; i++)
r[i].insert(i);
s  = f.RevSparseJac(m, r);
check_set.clear();
ok    &= s[0] == check_set;
check_set.insert(0);
check_set.insert(1);
ok    &= s[1] == check_set;
ok    &= s[2] == check_set;
// -------------------------------------------------------------------
// check forward Hessian sparsity for f_2 (x)
CPPAD_TESTVECTOR( std::set<size_t> ) r2(1), s2(1), h(n);
for(size_t j = 0; j < n; j++)
r2[0].insert(j);
s2[0].clear();
s2[0].insert(2);
h = f.ForSparseHes(r2, s2);
check_set.clear();
check_set.insert(0);
ok &= h[0] == check_set;
check_set.clear();
check_set.insert(1);
ok &= h[1] == check_set;
// -------------------------------------------------------------------
// check reverse Hessian sparsity for f_2 (x)
s3[0].clear();
s3[0].insert(2);
h = f.RevSparseHes(n, s3);
check_set.clear();
check_set.insert(0);
ok &= h[0] == check_set;
check_set.clear();
check_set.insert(1);
ok &= h[1] == check_set;
// -------------------------------------------------------------------
return ok;
}

Input File: example/atomic/eigen_mat_mul.cpp