1 # ifndef CPPAD_CORE_CHKPOINT_ONE_FOR_SPARSE_JAC_HPP
2 # define CPPAD_CORE_CHKPOINT_ONE_FOR_SPARSE_JAC_HPP
3 /* --------------------------------------------------------------------------
4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-19 Bradley M. Bell
5
6 CppAD is distributed under the terms of the
7 Eclipse Public License Version 2.0.
8
9 This Source Code may also be made available under the following
10 Secondary License when the conditions for such availability set forth
11 in the Eclipse Public License, Version 2.0 are satisfied:
12 GNU General Public License, Version 2.0 or later.
13 ---------------------------------------------------------------------------- */
14
15 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
16
17 template <class Base>
18 template <class sparsity_type>
for_sparse_jac(size_t q,const sparsity_type & r,sparsity_type & s,const vector<Base> & x)19 bool checkpoint<Base>::for_sparse_jac(
20 size_t q ,
21 const sparsity_type& r ,
22 sparsity_type& s ,
23 const vector<Base>& x )
24 { // make sure member_ is allocated for this thread
25 size_t thread = thread_alloc::thread_num();
26 allocate_member(thread);
27 //
28 // during user sparsity calculations
29 size_t m = member_[thread]->f_.Range();
30 size_t n = member_[thread]->f_.Domain();
31 if( member_[thread]->jac_sparse_bool_.size() == 0 )
32 set_jac_sparse_bool();
33 if( member_[thread]->jac_sparse_set_.n_set() != 0 )
34 member_[thread]->jac_sparse_set_.resize(0, 0);
35 CPPAD_ASSERT_UNKNOWN( member_[thread]->jac_sparse_bool_.size() == m * n );
36 CPPAD_ASSERT_UNKNOWN( member_[thread]->jac_sparse_set_.n_set() == 0 );
37 CPPAD_ASSERT_UNKNOWN( r.size() == n * q );
38 CPPAD_ASSERT_UNKNOWN( s.size() == m * q );
39 //
40 bool ok = true;
41 for(size_t i = 0; i < m; i++)
42 { for(size_t k = 0; k < q; k++)
43 s[i * q + k] = false;
44 }
45 // sparsity for s = jac_sparse_bool_ * r
46 for(size_t i = 0; i < m; i++)
47 { for(size_t k = 0; k < q; k++)
48 { // initialize sparsity for S(i,k)
49 bool s_ik = false;
50 // S(i,k) = sum_j J(i,j) * R(j,k)
51 for(size_t j = 0; j < n; j++)
52 { bool J_ij = member_[thread]->jac_sparse_bool_[ i * n + j];
53 bool R_jk = r[j * q + k ];
54 s_ik |= ( J_ij & R_jk );
55 }
56 s[i * q + k] = s_ik;
57 }
58 }
59 return ok;
60 }
61 template <class Base>
for_sparse_jac(size_t q,const vectorBool & r,vectorBool & s,const vector<Base> & x)62 bool checkpoint<Base>::for_sparse_jac(
63 size_t q ,
64 const vectorBool& r ,
65 vectorBool& s ,
66 const vector<Base>& x )
67 { // make sure member_ is allocated for this thread
68 size_t thread = thread_alloc::thread_num();
69 allocate_member(thread);
70 //
71 return for_sparse_jac< vectorBool >(q, r, s, x);
72 }
73 template <class Base>
for_sparse_jac(size_t q,const vector<bool> & r,vector<bool> & s,const vector<Base> & x)74 bool checkpoint<Base>::for_sparse_jac(
75 size_t q ,
76 const vector<bool>& r ,
77 vector<bool>& s ,
78 const vector<Base>& x )
79 { // make sure member_ is allocated for this thread
80 size_t thread = thread_alloc::thread_num();
81 allocate_member(thread);
82 //
83 return for_sparse_jac< vector<bool> >(q, r, s, x);
84 }
85 template <class Base>
for_sparse_jac(size_t q,const vector<std::set<size_t>> & r,vector<std::set<size_t>> & s,const vector<Base> & x)86 bool checkpoint<Base>::for_sparse_jac(
87 size_t q ,
88 const vector< std::set<size_t> >& r ,
89 vector< std::set<size_t> >& s ,
90 const vector<Base>& x )
91 { // make sure member_ is allocated for this thread
92 size_t thread = thread_alloc::thread_num();
93 allocate_member(thread);
94 //
95 // during user sparsity calculations
96 size_t m = member_[thread]->f_.Range();
97 size_t n = member_[thread]->f_.Domain();
98 if( member_[thread]->jac_sparse_bool_.size() != 0 )
99 member_[thread]->jac_sparse_bool_.clear();
100 if( member_[thread]->jac_sparse_set_.n_set() == 0 )
101 set_jac_sparse_set();
102 CPPAD_ASSERT_UNKNOWN( member_[thread]->jac_sparse_bool_.size() == 0 );
103 CPPAD_ASSERT_UNKNOWN( member_[thread]->jac_sparse_set_.n_set() == m );
104 CPPAD_ASSERT_UNKNOWN( member_[thread]->jac_sparse_set_.end() == n );
105 CPPAD_ASSERT_UNKNOWN( r.size() == n );
106 CPPAD_ASSERT_UNKNOWN( s.size() == m );
107
108 bool ok = true;
109 for(size_t i = 0; i < m; i++)
110 s[i].clear();
111
112 // sparsity for s = jac_sparse_set_ * r
113 for(size_t i = 0; i < m; i++)
114 { // compute row i of the return pattern
115 local::sparse::list_setvec::const_iterator set_itr(
116 member_[thread]->jac_sparse_set_, i
117 );
118 size_t j = *set_itr;
119 while(j < n )
120 { std::set<size_t>::const_iterator itr_j;
121 const std::set<size_t>& r_j( r[j] );
122 for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
123 { size_t k = *itr_j;
124 CPPAD_ASSERT_UNKNOWN( k < q );
125 s[i].insert(k);
126 }
127 j = *(++set_itr);
128 }
129 }
130
131 return ok;
132 }
133
134 } // END_CPPAD_NAMESPACE
135 # endif
136