Prev
Next
Index->
contents
reference
index
search
external
Up->
CppAD
AD
ADValued
atomic
atomic_base
atomic_mat_mul.cpp
atomic_mat_mul.hpp
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
atomic_mat_mul.hpp
Headings->
See Also
Matrix Dimensions
Start Class Definition
Constructor
Left Operand Element Index
Right Operand Element Index
Result Element Index
Forward Matrix Multiply
Reverse Matrix Multiply
forward
reverse
for_sparse_jac
rev_sparse_jac
rev_sparse_hes
End Class Definition
@(@\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}} }@)@Matrix Multiply as an Atomic Operation
See Also
atomic_eigen_mat_mul.hpp
Matrix Dimensions
This example puts the matrix dimensions in the atomic function arguments,
instead of the constructor
, so that they can
be different for different calls to the atomic function.
These dimensions are:
nr_left
number of rows in the left matrix
n_middle
rows in the left matrix and columns in right
nc_right
number of columns in the right matrix
Start Class Definition
# include <cppad/cppad.hpp>
namespace { // Begin empty namespace
using CppAD:: vector;
//
using CppAD:: set_union;
//
// matrix result = left * right
class atomic_mat_mul : public CppAD:: atomic_base< double > {
Constructor
public :
// ---------------------------------------------------------------------
// constructor
atomic_mat_mul ( void ) : CppAD:: atomic_base< double >( "mat_mul" )
{ }
private :
Left Operand Element Index
Index in the Taylor coefficient matrix
tx
of a left matrix element.
size_t left (
size_t i , // left matrix row index
size_t j , // left matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in tx
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert ( i < nr_left );
assert ( j < n_middle );
return ( 3 + i * n_middle + j) * nk + k;
}
Right Operand Element Index
Index in the Taylor coefficient matrix
tx
of a right matrix element.
size_t right (
size_t i , // right matrix row index
size_t j , // right matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in tx
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert ( i < n_middle );
assert ( j < nc_right );
size_t offset = 3 + nr_left * n_middle;
return ( offset + i * nc_right + j) * nk + k;
}
Result Element Index
Index in the Taylor coefficient matrix
ty
of a result matrix element.
size_t result (
size_t i , // result matrix row index
size_t j , // result matrix column index
size_t k , // Taylor coeffocient order
size_t nk , // number of Taylor coefficients in ty
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{ assert ( i < nr_left );
assert ( j < nc_right );
return ( i * nc_right + j) * nk + k;
}
Forward Matrix Multiply
Forward mode multiply Taylor coefficients in
tx
and sum into
ty
(for one pair of left and right orders)
void forward_multiply (
size_t k_left , // order for left coefficients
size_t k_right , // order for right coefficients
const vector< double >& tx , // domain space Taylor coefficients
vector< double >& ty , // range space Taylor coefficients
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t nk = tx. size () / nx;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
assert ( nk == ty. size () / ny );
# endif
//
size_t k_result = k_left + k_right;
assert ( k_result < nk );
//
for ( size_t i = 0 ; i < nr_left; i++)
{ for ( size_t j = 0 ; j < nc_right; j++)
{ double sum = 0.0 ;
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k_left, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k_right, nk, nr_left, n_middle, nc_right
);
sum += tx[ i_left] * tx[ i_right];
}
size_t i_result = result (
i, j, k_result, nk, nr_left, n_middle, nc_right
);
ty[ i_result] += sum;
}
}
}
Reverse Matrix Multiply
Reverse mode partials of Taylor coefficients and sum into
px
(for one pair of left and right orders)
void reverse_multiply (
size_t k_left , // order for left coefficients
size_t k_right , // order for right coefficients
const vector< double >& tx , // domain space Taylor coefficients
const vector< double >& ty , // range space Taylor coefficients
vector< double >& px , // partials w.r.t. tx
const vector< double >& py , // partials w.r.t. ty
size_t nr_left , // rows in left matrix
size_t n_middle , // rows in left and columns in right
size_t nc_right ) // columns in right matrix
{
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t nk = tx. size () / nx;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
assert ( nk == ty. size () / ny );
# endif
assert ( tx. size () == px. size () );
assert ( ty. size () == py. size () );
//
size_t k_result = k_left + k_right;
assert ( k_result < nk );
//
for ( size_t i = 0 ; i < nr_left; i++)
{ for ( size_t j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k_result, nk, nr_left, n_middle, nc_right
);
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k_left, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k_right, nk, nr_left, n_middle, nc_right
);
// sum += tx[i_left] * tx[i_right];
px[ i_left] += tx[ i_right] * py[ i_result];
px[ i_right] += tx[ i_left] * py[ i_result];
}
}
}
return ;
}
forward
Routine called by CppAD during Forward
mode.
virtual bool forward (
size_t q ,
size_t p ,
const vector< bool >& vx ,
vector< bool >& vy ,
const vector< double >& tx ,
vector< double >& ty
)
{ size_t n_order = p + 1 ;
size_t nr_left = size_t ( tx[ 0 * n_order + 0 ] );
size_t n_middle = size_t ( tx[ 1 * n_order + 0 ] );
size_t nc_right = size_t ( tx[ 2 * n_order + 0 ] );
# ifndef NDEBUG
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert ( vx. size () == 0 || nx == vx. size () );
assert ( vx. size () == 0 || ny == vy. size () );
assert ( nx * n_order == tx. size () );
assert ( ny * n_order == ty. size () );
size_t i, j, ell;
// check if we are computing vy information
if ( vx. size () > 0 )
{ size_t nk = 1 ;
size_t k = 0 ;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ bool var = false ;
for ( ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
bool nz_left = vx[ i_left] |( tx[ i_left] != 0 .);
bool nz_right = vx[ i_right]|( tx[ i_right] != 0 .);
// if not multiplying by the constant zero
if ( nz_left & nz_right )
var |= bool ( vx[ i_left]) | bool ( vx[ i_right]);
}
size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
vy[ i_result] = var;
}
}
}
// initialize result as zero
size_t k;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ for ( k = q; k <= p; k++)
{ size_t i_result = result (
i, j, k, n_order, nr_left, n_middle, nc_right
);
ty[ i_result] = 0.0 ;
}
}
}
for ( k = q; k <= p; k++)
{ // sum the produces that result in order k
for ( ell = 0 ; ell <= k; ell++)
forward_multiply (
ell, k - ell, tx, ty, nr_left, n_middle, nc_right
);
}
// all orders are implented, so always return true
return true ;
}
reverse
Routine called by CppAD during Reverse
mode.
virtual bool reverse (
size_t p ,
const vector< double >& tx ,
const vector< double >& ty ,
vector< double >& px ,
const vector< double >& py
)
{ size_t n_order = p + 1 ;
size_t nr_left = size_t ( tx[ 0 * n_order + 0 ] );
size_t n_middle = size_t ( tx[ 1 * n_order + 0 ] );
size_t nc_right = size_t ( tx[ 2 * n_order + 0 ] );
# ifndef NDEBUG
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert ( nx * n_order == tx. size () );
assert ( ny * n_order == ty. size () );
assert ( px. size () == tx. size () );
assert ( py. size () == ty. size () );
// initialize summation
for ( size_t i = 0 ; i < px. size (); i++)
px[ i] = 0.0 ;
// number of orders to differentiate
size_t k = n_order;
while ( k--)
{ // differentiate the produces that result in order k
for ( size_t ell = 0 ; ell <= k; ell++)
reverse_multiply (
ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right
);
}
// all orders are implented, so always return true
return true ;
}
for_sparse_jac
Routines called by CppAD during ForSparseJac
.
// boolean sparsity patterns
virtual bool for_sparse_jac (
size_t q ,
const vector< bool >& r ,
vector< bool >& s ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
# ifndef NDEBUG
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert ( nx == x. size () );
assert ( nx * q == r. size () );
assert ( ny * q == s. size () );
size_t p;
// sparsity for S(x) = f'(x) * R
size_t nk = 1 ;
size_t k = 0 ;
for ( size_t i = 0 ; i < nr_left; i++)
{ for ( size_t j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
for ( p = 0 ; p < q; p++)
s[ i_result * q + p] = false ;
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
for ( p = 0 ; p < q; p++)
{ // cast avoids Microsoft warning (should not be needed)
s[ i_result * q + p] |= bool ( r[ i_left * q + p ] );
s[ i_result * q + p] |= bool ( r[ i_right * q + p ] );
}
}
}
}
return true ;
}
// set sparsity patterns
virtual bool for_sparse_jac (
size_t q ,
const vector< std:: set< size_t> >& r ,
vector< std:: set< size_t> >& s ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
# ifndef NDEBUG
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
size_t ny = nr_left * nc_right;
# endif
assert ( nx == x. size () );
assert ( nx == r. size () );
assert ( ny == s. size () );
// sparsity for S(x) = f'(x) * R
size_t nk = 1 ;
size_t k = 0 ;
for ( size_t i = 0 ; i < nr_left; i++)
{ for ( size_t j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
s[ i_result]. clear ();
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
//
s[ i_result] = set_union ( s[ i_result], r[ i_left] );
s[ i_result] = set_union ( s[ i_result], r[ i_right] );
}
}
}
return true ;
}
rev_sparse_jac
Routines called by CppAD during RevSparseJac
.
// boolean sparsity patterns
virtual bool rev_sparse_jac (
size_t q ,
const vector< bool >& rt ,
vector< bool >& st ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
assert ( nx == x. size () );
assert ( nx * q == st. size () );
assert ( ny * q == rt. size () );
size_t i, j, p;
// initialize
for ( i = 0 ; i < nx; i++)
{ for ( p = 0 ; p < q; p++)
st[ i * q + p ] = false ;
}
// sparsity for S(x)^T = f'(x)^T * R^T
size_t nk = 1 ;
size_t k = 0 ;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
for ( p = 0 ; p < q; p++)
{ st[ i_left * q + p] |= bool ( rt[ i_result * q + p] );
st[ i_right* q + p] |= bool ( rt[ i_result * q + p] );
}
}
}
}
return true ;
}
// set sparsity patterns
virtual bool rev_sparse_jac (
size_t q ,
const vector< std:: set< size_t> >& rt ,
vector< std:: set< size_t> >& st ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
assert ( nx == x. size () );
assert ( nx == st. size () );
assert ( ny == rt. size () );
size_t i, j;
// initialize
for ( i = 0 ; i < nx; i++)
st[ i]. clear ();
// sparsity for S(x)^T = f'(x)^T * R^T
size_t nk = 1 ;
size_t k = 0 ;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
//
st[ i_left] = set_union ( st[ i_left], rt[ i_result]);
st[ i_right] = set_union ( st[ i_right], rt[ i_result]);
}
}
}
return true ;
}
rev_sparse_hes
Routines called by RevSparseHes
.
// set sparsity patterns
virtual bool rev_sparse_hes (
const vector< bool >& vx,
const vector< bool >& s ,
vector< bool >& t ,
size_t q ,
const vector< std:: set< size_t> >& r ,
const vector< std:: set< size_t> >& u ,
vector< std:: set< size_t> >& v ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
assert ( x. size () == nx );
assert ( vx. size () == nx );
assert ( t. size () == nx );
assert ( r. size () == nx );
assert ( v. size () == nx );
assert ( s. size () == ny );
assert ( u. size () == ny );
//
size_t i, j;
//
// initilaize sparsity patterns as false
for ( j = 0 ; j < nx; j++)
{ t[ j] = false ;
v[ j]. clear ();
}
size_t nk = 1 ;
size_t k = 0 ;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
//
// Compute sparsity for T(x) = S(x) * f'(x).
// We need not use vx with f'(x) back propagation.
t[ i_left] |= bool ( s[ i_result] );
t[ i_right] |= bool ( s[ i_result] );
// V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
// U(x) = g''(y) * f'(x) * R
// S(x) = g'(y)
// back propagate f'(x)^T * U(x)
// (no need to use vx with f'(x) propogation)
v[ i_left] = set_union ( v[ i_left], u[ i_result] );
v[ i_right] = set_union ( v[ i_right], u[ i_result] );
// back propagate S(x) * f''(x) * R
// (here is where we must check for cross terms)
if ( s[ i_result] & vx[ i_left] & vx[ i_right] )
{ v[ i_left] = set_union ( v[ i_left], r[ i_right] );
v[ i_right] = set_union ( v[ i_right], r[ i_left] );
}
}
}
}
return true ;
}
// bool sparsity
virtual bool rev_sparse_hes (
const vector< bool >& vx,
const vector< bool >& s ,
vector< bool >& t ,
size_t q ,
const vector< bool >& r ,
const vector< bool >& u ,
vector< bool >& v ,
const vector< double >& x )
{
size_t nr_left = size_t ( CppAD:: Integer ( x[ 0 ] ) );
size_t n_middle = size_t ( CppAD:: Integer ( x[ 1 ] ) );
size_t nc_right = size_t ( CppAD:: Integer ( x[ 2 ] ) );
size_t nx = 3 + ( nr_left + nc_right) * n_middle;
# ifndef NDEBUG
size_t ny = nr_left * nc_right;
# endif
assert ( x. size () == nx );
assert ( vx. size () == nx );
assert ( t. size () == nx );
assert ( r. size () == nx * q );
assert ( v. size () == nx * q );
assert ( s. size () == ny );
assert ( u. size () == ny * q );
size_t i, j, p;
//
// initilaize sparsity patterns as false
for ( j = 0 ; j < nx; j++)
{ t[ j] = false ;
for ( p = 0 ; p < q; p++)
v[ j * q + p] = false ;
}
size_t nk = 1 ;
size_t k = 0 ;
for ( i = 0 ; i < nr_left; i++)
{ for ( j = 0 ; j < nc_right; j++)
{ size_t i_result = result (
i, j, k, nk, nr_left, n_middle, nc_right
);
for ( size_t ell = 0 ; ell < n_middle; ell++)
{ size_t i_left = left (
i, ell, k, nk, nr_left, n_middle, nc_right
);
size_t i_right = right (
ell, j, k, nk, nr_left, n_middle, nc_right
);
//
// Compute sparsity for T(x) = S(x) * f'(x).
// We so not need to use vx with f'(x) propagation.
t[ i_left] |= bool ( s[ i_result] );
t[ i_right] |= bool ( s[ i_result] );
// V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
// U(x) = g''(y) * f'(x) * R
// S(x) = g'(y)
// back propagate f'(x)^T * U(x)
// (no need to use vx with f'(x) propogation)
for ( p = 0 ; p < q; p++)
{ v[ i_left * q + p] |= bool ( u[ i_result * q + p] );
v[ i_right * q + p] |= bool ( u[ i_result * q + p] );
}
// back propagate S(x) * f''(x) * R
// (here is where we must check for cross terms)
if ( s[ i_result] & vx[ i_left] & vx[ i_right] )
{ for ( p = 0 ; p < q; p++)
{ v[ i_left * q + p] |= bool ( r[ i_right * q + p] );
v[ i_right * q + p] |= bool ( r[ i_left * q + p] );
}
}
}
}
}
return true ;
}
End Class Definition
} ; // End of mat_mul class
} // End empty namespace
Input File: cppad/example/mat_mul.hpp