type_x
This example also demonstrates the use of the
type_x
argument to the atomic_three
virtual functions.
Matrix Dimensions
The first three components of the argument vector
ax
in the call
afun(ax, ay)
are parameters and contain the matrix dimensions.
This enables them to be different for each use of the same atomic
function
afun
.
These dimensions are:
ax[0]
nr_left
number of rows in the left matrix and result matrix
ax[1]
n_middle
columns in the left matrix and rows in right matrix
ax[2]
nc_right
number of columns in the right matrix and result matrix
Left Matrix
The number of elements in the left matrix is
n_left = nr_left * n_middle
The elements are in
ax[3]
through
ax[2+n_left]
in row major order.
Right Matrix
The number of elements in the right matrix is
n_right = n_middle * nc_right
The elements are in
ax[3+n_left]
through
ax[2+n_left+n_right]
in row major order.
Result Matrix
The number of elements in the result matrix is
n_result = nr_left * nc_right
The elements are in
ay[0]
through
ay[n_result-1]
in row major order.
# include <cppad/cppad.hpp>
namespace { // Begin empty namespaceusing CppAD::vector;
//// matrix result = left * rightclass atomic_mat_mul : public CppAD::atomic_three<double> {
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;
}
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 coefficientsconst 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 coefficientsconst vector<double>& tx , // domain space Taylor coefficientsconst vector<double>& ty , // range space Taylor coefficients
vector<double>& px , // partials w.r.t. txconst 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 );
# endifassert( 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;
}
