\(\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}} }\)
cppad_sparse_jacobian.cpp¶
View page sourceCppad Speed: Sparse Jacobian¶
Specifications¶
See link_sparse_jacobian .
Implementation¶
# include <cppad/cppad.hpp>
# include <cppad/speed/uniform_01.hpp>
# include <cppad/speed/sparse_jac_fun.hpp>
// Note that CppAD uses global_option["memory"] at the main program level
# include <map>
extern std::map<std::string, bool> global_option;
// see comments in main program for this external
extern size_t global_cppad_thread_alloc_inuse;
namespace {
using CppAD::vector;
typedef CppAD::AD<double> a_double;
typedef vector<size_t> s_vector;
typedef vector<bool> b_vector;
typedef vector<double> d_vector;
typedef vector<a_double> a_vector;
typedef CppAD::sparse_rc<s_vector> sparsity;
typedef CppAD::sparse_rcv<s_vector, d_vector> sparse_matrix;
void calc_sparsity(
CppAD::sparse_rc<s_vector>& pattern ,
CppAD::ADFun<double>& f )
{ bool reverse = global_option["revsparsity"];
bool transpose = false;
bool internal_bool = global_option["boolsparsity"];
bool dependency = false;
bool subgraph = global_option["subsparsity"];
size_t n = f.Domain();
size_t m = f.Range();
if( subgraph )
{ b_vector select_domain(n), select_range(m);
for(size_t j = 0; j < n; ++j)
select_domain[j] = true;
for(size_t i = 0; i < m; ++i)
select_range[i] = true;
f.subgraph_sparsity(
select_domain, select_range, transpose, pattern
);
}
else
{ size_t q = n;
if( reverse )
q = m;
//
CppAD::sparse_rc<s_vector> identity;
identity.resize(q, q, q);
for(size_t k = 0; k < q; k++)
identity.set(k, k, k);
//
if( reverse )
{ f.rev_jac_sparsity(
identity, transpose, dependency, internal_bool, pattern
);
}
else
{ f.for_jac_sparsity(
identity, transpose, dependency, internal_bool, pattern
);
}
}
}
// --------------------------------------------------------------------
void setup(
// inputs
size_t size ,
size_t m ,
const s_vector& row ,
const s_vector& col ,
// outputs
size_t& n_color ,
CppAD::ADFun<double>& f ,
sparse_matrix& subset ,
CppAD::sparse_jac_work& work )
{ // optimization options
std::string optimize_options =
"no_conditional_skip no_compare_op no_print_for_op";
if( global_option["val_graph"] )
optimize_options += " val_graph";
//
// default value for n_color
n_color = 0;
//
// independent variable vector
size_t nc = size;
a_vector a_x(nc);
d_vector x(nc);
//
// dependent variable vector
size_t nr = m;
a_vector a_y(nr);
//
// choose a value for independent variable vector
CppAD::uniform_01(nc, x);
for(size_t j = 0; j < nc; j++)
a_x[j] = x[j];
//
// declare independent variables
size_t abort_op_index = 0;
bool record_compare = false;
CppAD::Independent(a_x, abort_op_index, record_compare);
//
// AD computation of f(x)
size_t order = 0;
CppAD::sparse_jac_fun<a_double>(nr, nc, a_x, row, col, order, a_y);
//
// create function object f : x -> y
f.Dependent(a_x, a_y);
//
if( global_option["optimize"] )
f.optimize(optimize_options);
//
// coloring method
std::string coloring = "cppad";
# if CPPAD_HAS_COLPACK
if( global_option["colpack"] )
coloring = "colpack";
# else
CPPAD_ASSERT_UNKNOWN( ! global_option["colpack"] );
# endif
//
// sparsity pattern for subset of Jacobian that is evaluated
size_t nnz = row.size();
sparsity subset_pattern(nr, nc, nnz);
for(size_t k = 0; k < nnz; ++k)
subset_pattern.set(k, row[k], col[k]);
//
// sparse matrix for subset of Jacobian that is evaluated
subset = sparse_matrix( subset_pattern );
//
// maximum number of colors at once
size_t group_max = 25;
//
if( global_option["subgraph"] )
{ // This would cache some information in f, but would it enough ?
// The time it takes to compute derivatives that are not used
// slows down the test when onetape is false.
// f.subgraph_jac_rev(x, ac_subset);
}
else
{ // need full sparsity pattern
// (could use subset_sparsity, but pretend we do not konw that)
sparsity pattern;
calc_sparsity(pattern, f);
//
// Use forward mode to compute the Jacobian
// (this caches information in work),
work.clear();
n_color = f.sparse_jac_for(
group_max, x, subset, pattern, coloring, work
);
}
}
}
bool link_sparse_jacobian(
const std::string& job ,
size_t size ,
size_t repeat ,
size_t m ,
const CppAD::vector<size_t>& row ,
const CppAD::vector<size_t>& col ,
CppAD::vector<double>& x ,
CppAD::vector<double>& jacobian ,
size_t& n_color )
{ global_cppad_thread_alloc_inuse = 0;
// --------------------------------------------------------------------
// check global options
const char* valid[] = {
"memory", "onetape", "optimize", "subgraph",
"boolsparsity", "revsparsity", "subsparsity", "val_graph"
# if CPPAD_HAS_COLPACK
, "colpack"
# endif
};
size_t n_valid = sizeof(valid) / sizeof(valid[0]);
typedef std::map<std::string, bool>::iterator iterator;
//
for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
{ if( itr->second )
{ bool ok = false;
for(size_t i = 0; i < n_valid; i++)
ok |= itr->first == valid[i];
if( ! ok )
return false;
}
}
if( global_option["subsparsity"] )
{ if( global_option["boolsparisty"]
|| global_option["revsparsity"]
|| global_option["colpack"] )
return false;
}
// -----------------------------------------------------
// size corresponding to static_f
static size_t static_size = 0;
//
// function object corresponding to f(x)
static CppAD::ADFun<double> static_f;
//
// subset of Jacobian that we are using
static sparse_matrix static_subset;
//
// information used by for_sparse_jac_for
static CppAD::sparse_jac_work static_work;
//
// sparsity pattern not used because work is non-empty
sparsity empty_pattern;
// -----------------------------------------------------------------------
//
// default value for n_color
n_color = 0;
//
bool onetape = global_option["onetape"];
//
if( job == "setup" )
{ if( onetape )
{ setup(size, m, row, col,
n_color, static_f, static_subset, static_work
);
static_size = size;
}
else
{ static_size = 0;
}
return true;
}
if( job == "teardown" )
{ static_f = CppAD::ADFun<double>();
sparse_matrix empty_matrix;
static_subset.swap( empty_matrix );
static_work.clear();
static_size = 0;
return true;
}
// ------------------------------------------------------------------------
CPPAD_ASSERT_UNKNOWN( job == "run" );
//
// number of independent variables
static size_t n = size;
//
// maximum number of colors at once
size_t group_max = 25;
//
// coloring method
std::string coloring = "cppad";
if( global_option["colpack"] )
coloring = "colpack";
// ------------------------------------------------------
while(repeat--)
{ if( onetape )
{ if( size != static_size )
CPPAD_ASSERT_UNKNOWN( size == static_size );
}
else
{ setup(size, m, row, col,
n_color, static_f, static_subset, static_work
);
}
// choose a value for x
CppAD::uniform_01(n, x);
if( global_option["subgraph"] )
{ // user reverse mode because forward not yet implemented
static_f.subgraph_jac_rev(x, static_subset);
}
else
{ // Use forward mode because m > n (is this sufficient reason ?)
n_color = static_f.sparse_jac_for(group_max, x,
static_subset, empty_pattern, coloring, static_work
);
}
jacobian = static_subset.val();
}
size_t thread = CppAD::thread_alloc::thread_num();
global_cppad_thread_alloc_inuse = CppAD::thread_alloc::inuse(thread);
return true;
}