Prev Next Index-> contents reference index search external Up-> CppAD ADFun sparse_derivative sparse_jac sparse_jac_rev.cpp ADFun-> record_adfun drivers Forward Reverse sparsity_pattern sparse_derivative optimize abs_normal FunCheck check_for_nan sparse_derivative-> sparse_jac sparse_jacobian sparse_hes sparse_hessian subgraph_jac_rev sparse_jac-> sparse_jac_for.cpp sparse_jac_rev.cpp sparse_jac_rev.cpp Headings

$\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}} }$
Computing Sparse Jacobian Using Reverse Mode: Example and Test
# include <cppad/cppad.hpp> bool sparse_jac_rev(void) { bool ok = true; // using CppAD::AD; using CppAD::NearEqual; using CppAD::sparse_rc; using CppAD::sparse_rcv; // typedef CPPAD_TESTVECTOR(AD<double>) a_vector; typedef CPPAD_TESTVECTOR(double) d_vector; typedef CPPAD_TESTVECTOR(size_t) s_vector; // // domain space vector size_t n = 4; a_vector a_x(n); for(size_t j = 0; j < n; j++) a_x[j] = AD<double> (0); // // declare independent variables and starting recording CppAD::Independent(a_x); // size_t m = 3; a_vector a_y(m); a_y[0] = a_x[0] + a_x[1]; a_y[1] = a_x[2] + a_x[3]; a_y[2] = a_x[0] + a_x[1] + a_x[2] + a_x[3] * a_x[3] / 2.; // // create f: x -> y and stop tape recording CppAD::ADFun<double> f(a_x, a_y); // // new value for the independent variable vector d_vector x(n); for(size_t j = 0; j < n; j++) x[j] = double(j); /* [ 1 1 0 0 ] J(x) = [ 0 0 1 1 ] [ 1 1 1 x_3] */ // // row-major order values of J(x) size_t nnz = 8; s_vector check_row(nnz), check_col(nnz); d_vector check_val(nnz); for(size_t k = 0; k < nnz; k++) { // check_val if( k < 7 ) check_val[k] = 1.0; else check_val[k] = x[3]; // // check_row and check_col check_col[k] = k; if( k < 2 ) check_row[k] = 0; else if( k < 4 ) check_row[k] = 1; else { check_row[k] = 2; check_col[k] = k - 4; } } // // m by m identity matrix sparsity sparse_rc<s_vector> pattern_in(m, m, m); for(size_t k = 0; k < m; k++) pattern_in.set(k, k, k); // // sparsity for J(x) bool transpose = false; bool dependency = false; bool internal_bool = true; sparse_rc<s_vector> pattern_jac; f.rev_jac_sparsity( pattern_in, transpose, dependency, internal_bool, pattern_jac ); // // compute entire reverse mode Jacobian sparse_rcv<s_vector, d_vector> subset( pattern_jac ); CppAD::sparse_jac_work work; std::string coloring = "cppad"; size_t n_sweep = f.sparse_jac_rev(x, subset, pattern_jac, coloring, work); ok &= n_sweep == 2; // const s_vector row( subset.row() ); const s_vector col( subset.col() ); const d_vector val( subset.val() ); s_vector row_major = subset.row_major(); ok &= subset.nnz() == nnz; for(size_t k = 0; k < nnz; k++) { ok &= row[ row_major[k] ] == check_row[k]; ok &= col[ row_major[k] ] == check_col[k]; ok &= val[ row_major[k] ] == check_val[k]; } // // test using work stored by previous sparse_jac_rev sparse_rc<s_vector> pattern_not_used; std::string coloring_not_used; n_sweep = f.sparse_jac_rev(x, subset, pattern_jac, coloring, work); ok &= n_sweep == 2; for(size_t k = 0; k < nnz; k++) { ok &= row[ row_major[k] ] == check_row[k]; ok &= col[ row_major[k] ] == check_col[k]; ok &= val[ row_major[k] ] == check_val[k]; } // // compute non-zero in col 3 only, nr = m, nc = n, nnz = 2 sparse_rc<s_vector> pattern_col3(m, n, 2); pattern_col3.set(0, 1, 3); // row[0] = 1, col[0] = 3 pattern_col3.set(1, 2, 3); // row[1] = 2, col[1] = 3 sparse_rcv<s_vector, d_vector> subset_col3( pattern_col3 ); work.clear(); n_sweep = f.sparse_jac_rev(x, subset_col3, pattern_jac, coloring, work); ok &= n_sweep == 2; // const s_vector row_col3( subset_col3.row() ); const s_vector col_col3( subset_col3.col() ); const d_vector val_col3( subset_col3.val() ); ok &= subset_col3.nnz() == 2; // ok &= row_col3[0] == 1; ok &= col_col3[0] == 3; ok &= val_col3[0] == 1.0; // ok &= row_col3[1] == 2; ok &= col_col3[1] == 3; ok &= val_col3[1] == x[3]; // return ok; } 
Input File: example/sparse/sparse_jac_rev.cpp