\(\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_mat_mul_sparsity.cpp¶
View page sourceAtomic Matrix Multiply Sparsity Patterns: Example and Test¶
Purpose¶
This example demonstrates computing sparsity patterns with the atomic_four_mat_mul class.
f(x)¶
For a matrix \(A\) we define the function \(\R{rvec} ( A )\) to be the elements of \(A\) in row major order. For this example, the function \(f(x)\) is
Jacobian of f(x)¶
The Jacobian of \(f(x)\) is
Hessian¶
The function \(f_2 (x)\) is
The Hessian of \(f_2(x)\) is
where the first row is the column index, and the first column is the row index, for the corresponding matrix entries above.
Source¶
# include <cppad/cppad.hpp>
# include <cppad/example/atomic_four/mat_mul/mat_mul.hpp>
bool sparsity(void)
{ // ok, eps
bool ok = true;
//
// AD
using CppAD::AD;
using CppAD::sparse_rc;
// -----------------------------------------------------------------------
// Record f
// -----------------------------------------------------------------------
//
// afun
CppAD::atomic_mat_mul<double> afun("atomic_mat_mul");
//
// nleft, n_middle, n_right
size_t n_left = 2, n_middle = 2, n_right = 2;
//
// nx, ax
size_t nx = n_middle * (n_left + n_right);
CPPAD_TESTVECTOR( AD<double> ) ax(nx);
for(size_t j = 0; j < nx; ++j)
ax[j] = AD<double>(j + 2);
CppAD::Independent(ax);
//
// ny, ay
size_t ny = n_left * n_right;
CPPAD_TESTVECTOR( AD<double> ) ay(ny);
//
// ay
size_t call_id = afun.set(n_left, n_middle, n_right);
afun(call_id, ax, ay);
//
// f
CppAD::ADFun<double> f(ax, ay);
//
// s_vector
typedef CPPAD_TESTVECTOR(size_t) s_vector;
//
// eye_sparsity
// nx by nx identitty matrix
sparse_rc<s_vector> eye_sparsity;
eye_sparsity.resize(nx, nx, nx);
for(size_t i = 0; i < nx; ++i)
eye_sparsity.set(i, i, i);
//
// -----------------------------------------------------------------------
// jac_sparsity
bool transpose = false;
bool dependency = false;
bool internal_bool = false;
sparse_rc<s_vector> jac_sparsity;
f.for_jac_sparsity(
eye_sparsity, transpose, dependency, internal_bool, jac_sparsity
);
{ // check jac_sparsity
//
// row, col
const s_vector& row = jac_sparsity.row();
const s_vector& col = jac_sparsity.col();
s_vector row_major = jac_sparsity.row_major();
//
// ok
ok &= jac_sparsity.nnz() == 16;
for(size_t k = 0; k < jac_sparsity.nnz(); ++k)
ok &= row[ row_major[k] ] == k / 4;
// row 0
ok &= col[ row_major[0] ] == 0;
ok &= col[ row_major[1] ] == 1;
ok &= col[ row_major[2] ] == 4;
ok &= col[ row_major[3] ] == 6;
// row 1
ok &= col[ row_major[4] ] == 0;
ok &= col[ row_major[5] ] == 1;
ok &= col[ row_major[6] ] == 5;
ok &= col[ row_major[7] ] == 7;
// row 2
ok &= col[ row_major[8] ] == 2;
ok &= col[ row_major[9] ] == 3;
ok &= col[ row_major[10] ] == 4;
ok &= col[ row_major[11] ] == 6;
// row 3
ok &= col[ row_major[12] ] == 2;
ok &= col[ row_major[13] ] == 3;
ok &= col[ row_major[14] ] == 5;
ok &= col[ row_major[15] ] == 7;
}
// ----------------------------------------------------------------
//
// select_y
// corresponding to f_2
CPPAD_TESTVECTOR(bool) select_y(ny);
for(size_t i = 0; i < ny; ++i)
select_y[i] = false;
select_y[2] = true;
//
// hes_sparsity
transpose = false;
internal_bool = false;
sparse_rc<s_vector> hes_sparsity;
f.rev_hes_sparsity(select_y, transpose, internal_bool, hes_sparsity);
{ // check hes_sparsity
//
// row, col
const s_vector& row = hes_sparsity.row();
const s_vector& col = hes_sparsity.col();
s_vector row_major = hes_sparsity.row_major();
//
// ok
ok &= hes_sparsity.nnz() == 4;
//
ok &= row[ row_major[0] ] == 2;
ok &= col[ row_major[0] ] == 4;
//
ok &= row[ row_major[1] ] == 3;
ok &= col[ row_major[1] ] == 6;
//
ok &= row[ row_major[2] ] == 4;
ok &= col[ row_major[2] ] == 2;
//
ok &= row[ row_major[3] ] == 6;
ok &= col[ row_major[3] ] == 3;
}
return ok;
}