Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_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_mat_mul.cpp-> atomic_mat_mul.hpp Headings-> See Also Class Definition Use Atomic Function ---..Constructor ---..Recording ---..forward ---..reverse ---..option ---..for_sparse_jac ---..rev_sparse_jac ---..rev_sparse_hes

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

atomic_eigen_mat_mul.cpp

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

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

bool mat_mul(void)
{     bool ok = true;
size_t i, j;

Constructor

// -------------------------------------------------------------------
// object that multiplies  2 x 2  matrices
atomic_mat_mul afun;


Recording
     // start recording with four independent varables
size_t n = 4;
vector<double> x(n);
for(j = 0; j < n; j++)
ax[j] = x[j] = double(j + 1);

// ------------------------------------------------------------------
size_t nr_left = 2;
size_t n_middle  = 2;
size_t nc_right = 2;
vector< AD<double> > atom_x(3 + (nr_left + nc_right) * n_middle );

// matrix dimensions

// left matrix
atom_x[3] = ax[0];  // left[0, 0] = x0
atom_x[4] = ax[1];  // left[0, 1] = x1
atom_x[5] = 5.;     // left[1, 0] = 5
atom_x[6] = 6.;     // left[1, 1] = 6

// right matix
atom_x[7] = ax[2];  // right[0, 0] = x2
atom_x[8] = 7.;     // right[0, 1] = 7
atom_x[9] = ax[3];  // right[1, 0] = x3
atom_x[10] = 8.;     // right[1, 1] = 8
// ------------------------------------------------------------------
/*
[ x0 , x1 ] * [ x2 , 7 ] = [ x0*x2 + x1*x3 , x0*7 + x1*8 ]
[ 5  , 6  ]   [ x3 , 8 ]   [  5*x2 +  6*x3 ,  5*7 +  6*8 ]
*/
vector< AD<double> > atom_y(nr_left * nc_right);
afun(atom_x, atom_y);

ok &= (atom_y[0] == x[0]*x[2] + x[1]*x[3]) & Variable(atom_y[0]);
ok &= (atom_y[1] == x[0]*7.   + x[1]*8.  ) & Variable(atom_y[1]);
ok &= (atom_y[2] ==   5.*x[2] +   6.*x[3]) & Variable(atom_y[2]);
ok &= (atom_y[3] ==   5.*7.   +   6.*8.  ) & Parameter(atom_y[3]);

// ------------------------------------------------------------------
// define the function g : x -> atom_y
// g(x) = [ x0*x2 + x1*x3 , x0*7 + x1*8 , 5*x2  + 6*x3  , 5*7 + 6*8 ]^T
CppAD::ADFun<double> g(ax, atom_y);

forward
     // Test zero order forward mode evaluation of g(x)
size_t m = atom_y.size();
vector<double> y(m);
for(j = 0; j <  n; j++)
x[j] = double(j + 2);
y = g.Forward(0, x);
ok &= y[0] == x[0] * x[2] + x[1] * x[3];
ok &= y[1] == x[0] * 7.   + x[1] * 8.;
ok &= y[2] == 5. * x[2]   + 6. * x[3];
ok &= y[3] == 5. * 7.     + 6. * 8.;

//----------------------------------------------------------------------
// Test first order forward mode evaluation of g'(x) * [1, 2, 3, 4]^T
// g'(x) = [ x2, x3, x0, x1 ]
//         [ 7 ,  8,  0, 0  ]
//         [ 0 ,  0,  5, 6  ]
//         [ 0 ,  0,  0, 0  ]
for(j = 0; j <  n; j++)
dx[j] = double(j + 1);
dy = g.Forward(1, dx);
ok &= dy[0] == 1. * x[2] + 2. * x[3] + 3. * x[0] + 4. * x[1];
ok &= dy[1] == 1. * 7.   + 2. * 8.   + 3. * 0.   + 4. * 0.;
ok &= dy[2] == 1. * 0.   + 2. * 0.   + 3. * 5.   + 4. * 6.;
ok &= dy[3] == 1. * 0.   + 2. * 0.   + 3. * 0.   + 4. * 0.;

//----------------------------------------------------------------------
// Test second order forward mode
// g_0^2 (x) = [ 0, 0, 1, 0 ], g_0^2 (x) * [1] = [3]
//             [ 0, 0, 0, 1 ]              [2]   [4]
//             [ 1, 0, 0, 0 ]              [3]   [1]
//             [ 0, 1, 0, 0 ]              [4]   [2]
for(j = 0; j <  n; j++)
ddx[j] = 0.;
ddy = g.Forward(2, ddx);

// [1, 2, 3, 4] * g_0^2 (x) * [1, 2, 3, 4]^T = 1*3 + 2*4 + 3*1 + 4*2
ok &= 2. * ddy[0] == 1. * 3. + 2. * 4. + 3. * 1. + 4. * 2.;

// for i > 0, [1, 2, 3, 4] * g_i^2 (x) * [1, 2, 3, 4]^T = 0
ok &= ddy[1] == 0.;
ok &= ddy[2] == 0.;
ok &= ddy[3] == 0.;

reverse
     // Test second order reverse mode
for(i = 0; i < m; i++)
w[i] = 0.;
w[0] = 1.;
dw = g.Reverse(2, w);

// g_0'(x) = [ x2, x3, x0, x1 ]
ok &= dw[0*2 + 0] == x[2];
ok &= dw[1*2 + 0] == x[3];
ok &= dw[2*2 + 0] == x[0];
ok &= dw[3*2 + 0] == x[1];

// g_0'(x)   * [1, 2, 3, 4]  = 1 * x2 + 2 * x3 + 3 * x0 + 4 * x1
// g_0^2 (x) * [1, 2, 3, 4]  = [3, 4, 1, 2]
ok &= dw[0*2 + 1] == 3.;
ok &= dw[1*2 + 1] == 4.;
ok &= dw[2*2 + 1] == 1.;
ok &= dw[3*2 + 1] == 2.;

option
     //----------------------------------------------------------------------
// Test both the boolean and set sparsity at the atomic level
for(size_t sparse_index = 0; sparse_index < 2; sparse_index++)
{     if( sparse_index == 0 )
else     afun.option( CppAD::atomic_base<double>::set_sparsity_enum );

for_sparse_jac
     // Test forward Jacobian sparsity pattern
/*
g(x) = [ x0*x2 + x1*x3 , x0*7 + x1*8 , 5*x2  + 6*x3  , 5*7 + 6*8 ]^T
so the sparsity pattern should be
s[0] = {0, 1, 2, 3}
s[1] = {0, 1}
s[2] = {2, 3}
s[3] = {}
*/
for(j = 0; j <  n; j++)
{     assert( r[j].empty() );
r[j].insert(j);
}
s = g.ForSparseJac(n, r);
for(j = 0; j <  n; j++)
{     // s[0] = {0, 1, 2, 3}
ok &= s[0].find(j) != s[0].end();
// s[1] = {0, 1}
if( j == 0 || j == 1 )
ok &= s[1].find(j) != s[1].end();
else     ok &= s[1].find(j) == s[1].end();
// s[2] = {2, 3}
if( j == 2 || j == 3 )
ok &= s[2].find(j) != s[2].end();
else     ok &= s[2].find(j) == s[2].end();
}
// s[3] == {}
ok &= s[3].empty();

rev_sparse_jac
     // Test reverse Jacobian sparsity pattern
for(i = 0; i <  m; i++)
{     s[i].clear();
s[i].insert(i);
}
r = g.RevSparseJac(m, s);
for(j = 0; j <  n ; j++)
{     // r[0] = {0, 1, 2, 3}
ok &= r[0].find(j) != r[0].end();
// r[1] = {0, 1}
if( j == 0 || j == 1 )
ok &= r[1].find(j) != r[1].end();
else     ok &= r[1].find(j) == r[1].end();
// r[2] = {2, 3}
if( j == 2 || j == 3 )
ok &= r[2].find(j) != r[2].end();
else     ok &= r[2].find(j) == r[2].end();
}
// r[3] == {}
ok &= r[3].empty();

rev_sparse_hes
     /* Test reverse Hessian sparsity pattern
g_0^2 (x) = [ 0, 0, 1, 0 ] and for i > 0, g_i^2 = 0
[ 0, 0, 0, 1 ]
[ 1, 0, 0, 0 ]
[ 0, 1, 0, 0 ]
so for the sparsity pattern for the first component of g is
h[0] = {2}
h[1] = {3}
h[2] = {0}
h[3] = {1}
*/
t[0].clear();
t[0].insert(0);
h = g.RevSparseHes(n, t);
size_t check[] = {2, 3, 0, 1};
for(j = 0; j <  n; j++)
{     // h[j] = { check[j] }
for(i = 0; i < n; i++)
{     if( i == check[j] )
ok &= h[j].find(i) != h[j].end();
else     ok &= h[j].find(i) == h[j].end();
}
}
t[0].clear();
for( j = 1; j < n; j++)
t[0].insert(j);
h = g.RevSparseHes(n, t);
for(j = 0; j <  n; j++)
{     // h[j] = { }
for(i = 0; i < n; i++)
ok &= h[j].find(i) == h[j].end();
}

//-----------------------------------------------------------------
} // end for(size_t sparse_index  ...
//-----------------------------------------------------------------

return ok;
}

Input File: example/atomic/mat_mul.cpp