1 # ifndef CPPAD_CORE_SPARSE_JACOBIAN_HPP
2 # define CPPAD_CORE_SPARSE_JACOBIAN_HPP
18 # define CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION 64
364 template <
class VectorBase,
class VectorSet,
class VectorSize>
366 const VectorBase& x ,
367 VectorSet& p_transpose ,
368 const VectorSize& row ,
369 const VectorSize& col ,
386 CheckSimpleVector<Base, VectorBase>();
392 size_t K = size_t(jac.size());
403 if( color.
size() == 0 )
414 # if CPPAD_HAS_COLPACK
415 local::color_general_colpack(p_transpose, col, row, color);
419 "SparseJacobianForward: work.color_method = colpack "
420 "and colpack_prefix missing from cmake command line."
426 "SparseJacobianForward: work.color_method is not valid."
432 for(k = 0; k < K; k++)
433 key[k] = color[ col[k] ];
437 for(j = 0; j < n; j++)
if( color[j] < n )
438 n_color = std::max(n_color, color[j] + 1);
441 for(k = 0; k < K; k++)
444 # if CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION == 1
446 VectorBase dx(n), dy(m);
450 for(ell = 0; ell < n_color; ell++)
454 for(j = 0; j < n; j++)
456 if( color[j] == ell )
463 while( k < K && color[ col[order[k]] ] == ell )
464 { jac[ order[k] ] = dy[row[order[k]]];
474 size_t count_color = 0;
477 while( count_color < n_color )
479 size_t r = std::min(max_r , n_color - count_color);
480 VectorBase dx(n * r), dy(m * r);
483 for(ell = 0; ell < r; ell++)
485 for(j = 0; j < n; j++)
486 { dx[j * r + ell] = zero;
487 if( color[j] == ell + count_color )
488 dx[j * r + ell] = one;
492 dy = Forward(q, r, dx);
495 for(ell = 0; ell < r; ell++)
497 while( k < K && color[ col[order[k]] ] == ell + count_color )
498 { jac[ order[k] ] = dy[ row[order[k]] * r + ell ];
561 template <
class VectorBase,
class VectorSet,
class VectorSize>
563 const VectorBase& x ,
565 const VectorSize& row ,
566 const VectorSize& col ,
583 CheckSimpleVector<Base, VectorBase>();
589 size_t K = size_t(jac.size());
600 if( color.
size() == 0 )
611 # if CPPAD_HAS_COLPACK
612 local::color_general_colpack(p, row, col, color);
616 "SparseJacobianReverse: work.color_method = colpack "
617 "and colpack_prefix missing from cmake command line."
623 "SparseJacobianReverse: work.color_method is not valid."
629 for(k = 0; k < K; k++)
630 key[k] = color[ row[k] ];
634 for(i = 0; i < m; i++)
if( color[i] < m )
635 n_color = std::max(n_color, color[i] + 1);
644 for(k = 0; k < K; k++)
649 for(ell = 0; ell < n_color; ell++)
653 for(i = 0; i < m; i++)
655 if( color[i] == ell )
662 while( k < K && color[ row[order[k]] ] == ell )
663 { jac[ order[k] ] = dw[col[order[k]]];
727 template <
class VectorBase,
class VectorSet,
class VectorSize>
729 const VectorBase& x ,
731 const VectorSize& row ,
732 const VectorSize& col ,
738 size_t K = jac.size();
742 size_t(x.size()) == n ,
743 "SparseJacobianForward: size of x not equal domain dimension for f."
746 size_t(row.size()) == K &&
size_t(col.size()) == K ,
747 "SparseJacobianForward: either r or c does not have "
748 "the same size as jac."
752 "SparseJacobianForward: invalid value in work."
754 for(k = 0; k < K; k++)
757 "SparseJacobianForward: invalid value in r."
761 "SparseJacobianForward: invalid value in c."
767 "SparseJacobianForward: invalid value in work."
777 Pattern_type s_transpose;
779 {
bool transpose =
true;
780 const char* error_msg =
"SparseJacobianForward: transposed sparsity"
781 " pattern does not have proper row or column dimension";
784 n_sweep = SparseJacobianFor(x, s_transpose, row, col, jac, work);
842 template <
class VectorBase,
class VectorSet,
class VectorSize>
844 const VectorBase& x ,
846 const VectorSize& row ,
847 const VectorSize& col ,
853 size_t K = jac.size();
857 size_t(x.size()) == n ,
858 "SparseJacobianReverse: size of x not equal domain dimension for f."
861 size_t(row.size()) == K &&
size_t(col.size()) == K ,
862 "SparseJacobianReverse: either r or c does not have "
863 "the same size as jac."
867 "SparseJacobianReverse: invalid value in work."
869 for(k = 0; k < K; k++)
872 "SparseJacobianReverse: invalid value in r."
876 "SparseJacobianReverse: invalid value in c."
882 "SparseJacobianReverse: invalid value in work."
894 {
bool transpose =
false;
895 const char* error_msg =
"SparseJacobianReverse: sparsity"
896 " pattern does not have proper row or column dimension";
899 n_sweep = SparseJacobianRev(x, s, row, col, jac, work);
931 template <
class Base>
932 template <
class VectorBase,
class VectorSet>
934 const VectorBase& x,
const VectorSet& p
940 VectorBase jac(m * n);
943 size_t(x.size()) == n,
944 "SparseJacobian: size of x not equal domain size for f."
946 CheckSimpleVector<Base, VectorBase>();
953 for(i = 0; i < m; i++)
954 for(j = 0; j < n; j++)
955 jac[i * n + j] = zero;
963 Pattern_type s_transpose;
964 bool transpose =
true;
965 const char* error_msg =
"SparseJacobian: transposed sparsity"
966 " pattern does not have proper row or column dimension";
970 for(j = 0; j < n; j++)
971 {
typename Pattern_type::const_iterator itr(s_transpose, j);
973 while( i != s_transpose.end() )
984 SparseJacobianFor(x, s_transpose, row, col, J, work);
987 for(k = 0; k < K; k++)
988 jac[ row[k] * n + col[k] ] = J[k];
994 bool transpose =
false;
995 const char* error_msg =
"SparseJacobian: sparsity"
996 " pattern does not have proper row or column dimension";
1000 for(i = 0; i < m; i++)
1001 {
typename Pattern_type::const_iterator itr(s, i);
1003 while( j != s.end() )
1014 SparseJacobianRev(x, s, row, col, J, work);
1017 for(k = 0; k < K; k++)
1018 jac[ row[k] * n + col[k] ] = J[k];
1046 template <
class Base>
1047 template <
class VectorBase>
1052 size_t n = Domain();
1055 VectorBool p(m * n);
1061 VectorBool r(n * n);
1062 for(j = 0; j < n; j++)
1063 {
for(k = 0; k < n; k++)
1064 r[j * n + k] =
false;
1065 r[j * n + j] =
true;
1067 p = ForSparseJac(n, r);
1073 VectorBool s(m * m);
1074 for(i = 0; i < m; i++)
1075 {
for(k = 0; k < m; k++)
1076 s[i * m + k] =
false;
1077 s[i * m + i] =
true;
1079 p = RevSparseJac(m, s);
1081 return SparseJacobian(x, p);
1085 # undef CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
sparse_jacobian_work(void)
constructor
size_t SparseJacobianReverse(const VectorBase &x, const VectorSet &p, const VectorSize &r, const VectorSize &c, VectorBase &jac, sparse_jacobian_work &work)
Compute user specified subset of a sparse Jacobian using forward mode.
void color_general_cppad(const VectorSet &pattern, const VectorSize &row, const VectorSize &col, CppAD::vector< size_t > &color)
Determine which rows of a general sparse matrix can be computed together; i.e., do not have non-zero ...
std::string color_method
Coloring method: "cppad", or "colpack" (this field is set by user)
CppAD::vector< size_t > color
results of the coloring algorithm
void clear(void)
free memory and set number of elements to zero
VectorBase SparseJacobian(const VectorBase &x)
calculate sparse Jacobians
size_t SparseJacobianForward(const VectorBase &x, const VectorSet &p, const VectorSize &r, const VectorSize &c, VectorBase &jac, sparse_jacobian_work &work)
Compute user specified subset of a sparse Jacobian using forward mode.
Coloring algorithm for a general sparse matrix.
void resize(size_t n)
change the number of elements in this vector.
#define CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION
size_t size(void) const
number of elements currently in this vector.
void clear(void)
reset coloring method to its default and inform CppAD that color and order need to be recomputed ...
Two constant standard sets (currently used for concept checking).
void index_sort(const VectorKey &keys, VectorSize &ind)
Compute the indices that sort a vector of keys.
size_t SparseJacobianFor(const VectorBase &x, VectorSet &p_transpose, const VectorSize &row, const VectorSize &col, VectorBase &jac, sparse_jacobian_work &work)
Private helper function forward mode cases.
void push_back(const Type &s)
add an element to the back of this vector
class used by SparseJacobian to hold information so it does not need to be recomputed.
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
size_t SparseJacobianRev(const VectorBase &x, VectorSet &p, const VectorSize &row, const VectorSize &col, VectorBase &jac, sparse_jacobian_work &work)
Private helper function for reverse mode cases.
CppAD::vector< size_t > order
indices that sort the user row and col arrays by color
void sparsity_user2internal(sparse_list &internal, const VectorSet &user, size_t n_set, size_t end, bool transpose, const char *error_msg)
Copy a user vector of sets sparsity pattern to an internal sparse_list object.
Template structure used obtain the internal sparsity pattern type form the corresponding element type...