1 # ifndef CPPAD_CORE_FOR_JAC_SPARSITY_HPP
2 # define CPPAD_CORE_FOR_JAC_SPARSITY_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 $begin for_jac_sparsity$$
16 $spell
17 Jacobian
18 jac
19 bool
20 const
21 rc
22 cpp
23 $$
24
25 $section Forward Mode Jacobian Sparsity Patterns$$
26
27 $head Syntax$$
28 $icode%f%.for_jac_sparsity(
29 %pattern_in%, %transpose%, %dependency%, %internal_bool%, %pattern_out%
30 )%$$
31
32 $head Purpose$$
33 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
34 $cref/AD function/glossary/AD Function/$$ corresponding to
35 the operation sequence stored in $icode f$$.
36 Fix $latex R \in \B{R}^{n \times \ell}$$ and define the function
37 $latex \[
38 J(x) = F^{(1)} ( x ) * R
39 \] $$
40 Given the $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
41 $code for_jac_sparsity$$ computes a sparsity pattern for $latex J(x)$$.
42
43 $head x$$
44 Note that the sparsity pattern $latex J(x)$$ corresponds to the
45 operation sequence stored in $icode f$$ and does not depend on
46 the argument $icode x$$.
47 (The operation sequence may contain
48 $cref CondExp$$ and $cref VecAD$$ operations.)
49
50 $head SizeVector$$
51 The type $icode SizeVector$$ is a $cref SimpleVector$$ class with
52 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
53 $code size_t$$.
54
55 $head f$$
56 The object $icode f$$ has prototype
57 $codei%
58 ADFun<%Base%> %f%
59 %$$
60 The $cref ADFun$$ object $icode f$$ is not $code const$$.
61 After a call to $code for_jac_sparsity$$, a sparsity pattern
62 for each of the variables in the operation sequence
63 is held in $icode f$$ for possible later use during
64 reverse Hessian sparsity calculations.
65
66 $subhead size_forward_bool$$
67 After $code for_jac_sparsity$$, if $icode k$$ is a $code size_t$$ object,
68 $codei%
69 %k% = %f%.size_forward_bool()
70 %$$
71 sets $icode k$$ to the amount of memory (in unsigned character units)
72 used to store the
73 $cref/boolean vector/glossary/Sparsity Pattern/Boolean Vector/$$
74 sparsity patterns.
75 If $icode internal_bool$$ if false, $icode k$$ will be zero.
76 Otherwise it will be non-zero.
77 If you do not need this information for $cref RevSparseHes$$
78 calculations, it can be deleted
79 (and the corresponding memory freed) using
80 $codei%
81 %f%.size_forward_bool(0)
82 %$$
83 after which $icode%f%.size_forward_bool()%$$ will return zero.
84
85 $subhead size_forward_set$$
86 After $code for_jac_sparsity$$, if $icode k$$ is a $code size_t$$ object,
87 $codei%
88 %k% = %f%.size_forward_set()
89 %$$
90 sets $icode k$$ to the amount of memory (in unsigned character units)
91 used to store the
92 $cref/vector of sets/glossary/Sparsity Pattern/Vector of Sets/$$
93 sparsity patterns.
94 If $icode internal_bool$$ if true, $icode k$$ will be zero.
95 Otherwise it will be non-zero.
96 If you do not need this information for future $cref rev_hes_sparsity$$
97 calculations, it can be deleted
98 (and the corresponding memory freed) using
99 $codei%
100 %f%.size_forward_set(0)
101 %$$
102 after which $icode%f%.size_forward_set()%$$ will return zero.
103
104 $head pattern_in$$
105 The argument $icode pattern_in$$ has prototype
106 $codei%
107 const sparse_rc<%SizeVector%>& %pattern_in%
108 %$$
109 see $cref sparse_rc$$.
110 If $icode transpose$$ it is false (true),
111 $icode pattern_in$$ is a sparsity pattern for $latex R$$ ($latex R^\R{T}$$).
112
113 $head transpose$$
114 This argument has prototype
115 $codei%
116 bool %transpose%
117 %$$
118 See $cref/pattern_in/for_jac_sparsity/pattern_in/$$ above and
119 $cref/pattern_out/for_jac_sparsity/pattern_out/$$ below.
120
121 $head dependency$$
122 This argument has prototype
123 $codei%
124 bool %dependency%
125 %$$
126 see $cref/pattern_out/for_jac_sparsity/pattern_out/$$ below.
127
128 $head internal_bool$$
129 This argument has prototype
130 $codei%
131 bool %internal_bool%
132 %$$
133 If this is true, calculations are done with sets represented by a vector
134 of boolean values. Otherwise, a vector of sets of integers is used.
135
136 $head pattern_out$$
137 This argument has prototype
138 $codei%
139 sparse_rc<%SizeVector%>& %pattern_out%
140 %$$
141 This input value of $icode pattern_out$$ does not matter.
142 If $icode transpose$$ it is false (true),
143 upon return $icode pattern_out$$ is a sparsity pattern for
144 $latex J(x)$$ ($latex J(x)^\R{T}$$).
145 If $icode dependency$$ is true, $icode pattern_out$$ is a
146 $cref/dependency pattern/dependency.cpp/Dependency Pattern/$$
147 instead of sparsity pattern.
148
149 $head Sparsity for Entire Jacobian$$
150 Suppose that
151 $latex R$$ is the $latex n \times n$$ identity matrix.
152 In this case, $icode pattern_out$$ is a sparsity pattern for
153 $latex F^{(1)} ( x )$$ ( $latex F^{(1)} (x)^\R{T}$$ )
154 if $icode transpose$$ is false (true).
155
156 $head Example$$
157 $children%
158 example/sparse/for_jac_sparsity.cpp
159 %$$
160 The file
161 $cref for_jac_sparsity.cpp$$
162 contains an example and test of this operation.
163
164 $end
165 -----------------------------------------------------------------------------
166 */
167 # include <cppad/core/ad_fun.hpp>
168 # include <cppad/local/sparse/internal.hpp>
169
170 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
171
172 /*!
173 Forward Jacobian sparsity patterns.
174
175 \tparam Base
176 is the base type for this recording.
177
178 \tparam SizeVector
179 is the simple vector with elements of type size_t that is used for
180 row, column index sparsity patterns.
181
182 \param pattern_in
183 is the sparsity pattern for for R or R^T depending on transpose.
184
185 \param transpose
186 Is the input and returned sparsity pattern transposed.
187
188 \param dependency
189 Are the derivatives with respect to left and right of the expression below
190 considered to be non-zero:
191 \code
192 CondExpRel(left, right, if_true, if_false)
193 \endcode
194 This is used by the optimizer to obtain the correct dependency relations.
195
196 \param internal_bool
197 If this is true, calculations are done with sets represented by a vector
198 of boolean values. Othewise, a vector of standard sets is used.
199
200 \param pattern_out
201 The value of transpose is false (true),
202 the return value is a sparsity pattern for J(x) ( J(x)^T ) where
203 \f[
204 J(x) = F^{(1)} (x) * R
205 \f]
206 Here F is the function corresponding to the operation sequence
207 and x is any argument value.
208 */
209 template <class Base, class RecBase>
210 template <class SizeVector>
for_jac_sparsity(const sparse_rc<SizeVector> & pattern_in,bool transpose,bool dependency,bool internal_bool,sparse_rc<SizeVector> & pattern_out)211 void ADFun<Base,RecBase>::for_jac_sparsity(
212 const sparse_rc<SizeVector>& pattern_in ,
213 bool transpose ,
214 bool dependency ,
215 bool internal_bool ,
216 sparse_rc<SizeVector>& pattern_out )
217 {
218 // used to identify the RecBase type in calls to sweeps
219 RecBase not_used_rec_base(0.0);
220 //
221 // number or rows, columns, and non-zeros in pattern_in
222 size_t nr_in = pattern_in.nr();
223 size_t nc_in = pattern_in.nc();
224 //
225 size_t n = nr_in;
226 size_t ell = nc_in;
227 if( transpose )
228 std::swap(n, ell);
229 //
230 CPPAD_ASSERT_KNOWN(
231 n == Domain() ,
232 "for_jac_sparsity: number rows in R "
233 "is not equal number of independent variables."
234 );
235 bool zero_empty = true;
236 bool input_empty = true;
237 if( internal_bool )
238 { // allocate memory for bool sparsity calculation
239 // (sparsity pattern is emtpy after a resize)
240 for_jac_sparse_pack_.resize(num_var_tape_, ell);
241 for_jac_sparse_set_.resize(0, 0);
242 //
243 // set sparsity patttern for independent variables
244 local::sparse::set_internal_pattern(
245 zero_empty ,
246 input_empty ,
247 transpose ,
248 ind_taddr_ ,
249 for_jac_sparse_pack_ ,
250 pattern_in
251 );
252
253 // compute sparsity for other variables
254 local::sweep::for_jac<addr_t>(
255 &play_,
256 dependency,
257 n,
258 num_var_tape_,
259 for_jac_sparse_pack_,
260 not_used_rec_base
261
262 );
263 // set the output pattern
264 local::sparse::get_internal_pattern(
265 transpose, dep_taddr_, for_jac_sparse_pack_, pattern_out
266 );
267 }
268 else
269 {
270 // allocate memory for set sparsity calculation
271 // (sparsity pattern is emtpy after a resize)
272 for_jac_sparse_set_.resize(num_var_tape_, ell);
273 for_jac_sparse_pack_.resize(0, 0);
274 //
275 // set sparsity patttern for independent variables
276 local::sparse::set_internal_pattern(
277 zero_empty ,
278 input_empty ,
279 transpose ,
280 ind_taddr_ ,
281 for_jac_sparse_set_ ,
282 pattern_in
283 );
284
285 // compute sparsity for other variables
286 local::sweep::for_jac<addr_t>(
287 &play_,
288 dependency,
289 n,
290 num_var_tape_,
291 for_jac_sparse_set_,
292 not_used_rec_base
293
294 );
295 // get the ouput pattern
296 local::sparse::get_internal_pattern(
297 transpose, dep_taddr_, for_jac_sparse_set_, pattern_out
298 );
299 }
300 return;
301 }
302
303
304 } // END_CPPAD_NAMESPACE
305 # endif
306