CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_JACOBIAN_INCLUDED 00003 # define CPPAD_SPARSE_JACOBIAN_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 /* 00016 $begin sparse_jacobian$$ 00017 $spell 00018 cppad 00019 colpack 00020 cmake 00021 recomputed 00022 valarray 00023 std 00024 CppAD 00025 Bool 00026 jac 00027 Jacobian 00028 const 00029 Taylor 00030 $$ 00031 00032 $section Sparse Jacobian: Easy Driver$$ 00033 $index SparseJacobian$$ 00034 $index jacobian, sparse$$ 00035 00036 $head Syntax$$ 00037 $icode%jac% = %f%.SparseJacobian(%x%) 00038 %jac% = %f%.SparseJacobian(%x%, %p%) 00039 %n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%) 00040 %n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%) 00041 %$$ 00042 00043 $head Purpose$$ 00044 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size, 00045 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$. 00046 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the 00047 $cref/AD function/glossary/AD Function/$$ 00048 corresponding to $icode f$$. 00049 The syntax above sets $icode jac$$ to the Jacobian 00050 $latex \[ 00051 jac = F^{(1)} (x) 00052 \] $$ 00053 This routine takes advantage of the sparsity of the Jacobian 00054 in order to reduce the amount of computation necessary. 00055 If $icode row$$ and $icode col$$ are present, it also takes 00056 advantage of the reduced set of elements of the Jacobian that 00057 need to be computed. 00058 One can use speed tests (e.g. $cref speed_test$$) 00059 to verify that results are computed faster 00060 than when using the routine $cref Jacobian$$. 00061 00062 $head f$$ 00063 The object $icode f$$ has prototype 00064 $codei% 00065 ADFun<%Base%> %f% 00066 %$$ 00067 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$ 00068 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below). 00069 00070 $head x$$ 00071 The argument $icode x$$ has prototype 00072 $codei% 00073 const %VectorBase%& %x% 00074 %$$ 00075 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below) 00076 and its size 00077 must be equal to $icode n$$, the dimension of the 00078 $cref/domain/seq_property/Domain/$$ space for $icode f$$. 00079 It specifies 00080 that point at which to evaluate the Jacobian. 00081 00082 $head p$$ 00083 The argument $icode p$$ is optional and has prototype 00084 $codei% 00085 const %VectorSet%& %p% 00086 %$$ 00087 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below). 00088 If it has elements of type $code bool$$, 00089 its size is $latex m * n$$. 00090 If it has elements of type $code std::set<size_t>$$, 00091 its size is $latex m$$ and all its set elements are between 00092 zero and $latex n - 1$$. 00093 It specifies a 00094 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00095 for the Jacobian $latex F^{(1)} (x)$$. 00096 $pre 00097 00098 $$ 00099 If this sparsity pattern does not change between calls to 00100 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 00101 (using $cref ForSparseJac$$ or $cref RevSparseJac$$) 00102 and then pass $icode p$$ to $codei SparseJacobian$$. 00103 In addition, 00104 if you specify $icode p$$, CppAD will use the same 00105 type of sparsity representation 00106 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00107 for its internal calculations. 00108 Otherwise, the representation 00109 for the internal calculations is unspecified. 00110 00111 $head row, col$$ 00112 The arguments $icode row$$ and $icode col$$ are optional and have prototype 00113 $codei% 00114 const %VectorSize%& %row% 00115 const %VectorSize%& %col% 00116 %$$ 00117 (see $cref/VectorSize/sparse_jacobian/VectorSize/$$ below). 00118 They specify which rows and columns of $latex F^{(1)} (x)$$ are 00119 computes and in what order. 00120 Not all the non-zero entries in $latex F^{(1)} (x)$$ need be computed, 00121 but all the entries specified by $icode row$$ and $icode col$$ 00122 must be possibly non-zero in the sparsity pattern. 00123 We use $latex K$$ to denote the value $icode%jac%.size()%$$ 00124 which must also equal the size of $icode row$$ and $icode col$$. 00125 Furthermore, 00126 for $latex k = 0 , \ldots , K-1$$, it must hold that 00127 $latex row[k] < m$$ and $latex col[k] < n$$. 00128 00129 $head jac$$ 00130 The result $icode jac$$ has prototype 00131 $codei% 00132 %VectorBase%& %jac% 00133 %$$ 00134 In the case where the arguments $icode row$$ and $icode col$$ are not present, 00135 the size of $icode jac$$ is $latex m * n$$ and 00136 for $latex i = 0 , \ldots , m-1$$, 00137 $latex j = 0 , \ldots , n-1$$, 00138 $latex \[ 00139 jac [ i * n + j ] = \D{ F_i }{ x_j } (x) 00140 \] $$ 00141 $pre 00142 00143 $$ 00144 In the case where the arguments $icode row$$ and $icode col$$ are present, 00145 we use $latex K$$ to denote the size of $icode jac$$. 00146 The input value of its elements does not matter. 00147 Upon return, for $latex k = 0 , \ldots , K - 1$$, 00148 $latex \[ 00149 jac [ k ] = \D{ F_i }{ x_j } (x) 00150 \; , \; 00151 \; {\rm where} \; 00152 i = row[k] 00153 \; {\rm and } \; 00154 j = col[k] 00155 \] $$ 00156 00157 $head work$$ 00158 If this argument is present, it has prototype 00159 $codei% 00160 sparse_jacobian_work& %work% 00161 %$$ 00162 This object can only be used with the routines 00163 $code SparseJacobianForward$$ and $code SparseJacobianReverse$$. 00164 During its the first use, information is stored in $icode work$$. 00165 This is used to reduce the work done by future calls to the same mode 00166 (forward or reverse), 00167 the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$. 00168 If a future call is for a different mode, 00169 or any of these values have changed, 00170 you must first call $icode%work%.clear()%$$ 00171 to inform CppAD that this information needs to be recomputed. 00172 00173 $subhead color_method$$ 00174 The coloring algorithm determines which columns (forward mode) 00175 or rows (reverse mode) can be computed during the same sweep. 00176 This field has prototype 00177 $codep% 00178 std::string %work%.color_method 00179 %$$ 00180 and its default value (after a constructor or $code clear()$$) 00181 is $code "cppad"$$. 00182 If $cref colpack_prefix$$ is specified on the 00183 $cref/cmake command/cmake/CMake Command/$$ line, 00184 you can set this method to $code "colpack"$$. 00185 This value only matters on the first call to $code sparse_jacobian$$ 00186 after the $icode work$$ constructor or a call to $code clear$$. 00187 00188 $head n_sweep$$ 00189 The return value $icode n_sweep$$ has prototype 00190 $codei% 00191 size_t %n_sweep% 00192 %$$ 00193 If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used, 00194 $icode n_sweep$$ is the number of first order forward (reverse) sweeps 00195 used to compute the requested Jacobian values. 00196 (This is also the number of colors determined by the coloring method 00197 mentioned above). 00198 This is proportional to the total work that $code SparseJacobian$$ does, 00199 not counting the zero order forward sweep, 00200 or the work to combine multiple columns (rows) into a single sweep. 00201 00202 $head VectorBase$$ 00203 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with 00204 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00205 $icode Base$$. 00206 The routine $cref CheckSimpleVector$$ will generate an error message 00207 if this is not the case. 00208 00209 $head VectorSet$$ 00210 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with 00211 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00212 $code bool$$ or $code std::set<size_t>$$; 00213 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00214 of the difference. 00215 The routine $cref CheckSimpleVector$$ will generate an error message 00216 if this is not the case. 00217 00218 $subhead Restrictions$$ 00219 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00220 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00221 corresponding set. 00222 According to section 26.3.2.3 of the 1998 C++ standard, 00223 $code std::valarray< std::set<size_t> >$$ does not satisfy 00224 this condition. 00225 00226 $head VectorSize$$ 00227 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with 00228 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00229 $code size_t$$. 00230 The routine $cref CheckSimpleVector$$ will generate an error message 00231 if this is not the case. 00232 00233 $head Uses Forward$$ 00234 After each call to $cref Forward$$, 00235 the object $icode f$$ contains the corresponding 00236 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00237 After a call to any of the sparse Jacobian routines, 00238 the zero order Taylor coefficients correspond to 00239 $icode%f%.Forward(0, %x%)%$$ 00240 and the other coefficients are unspecified. 00241 00242 After $code SparseJacobian$$, 00243 the previous calls to $cref Forward$$ are undefined. 00244 00245 $head Example$$ 00246 $children% 00247 example/sparse_jacobian.cpp 00248 %$$ 00249 The routine 00250 $cref sparse_jacobian.cpp$$ 00251 is examples and tests of $code sparse_jacobian$$. 00252 It return $code true$$, if it succeeds and $code false$$ otherwise. 00253 00254 $end 00255 ============================================================================== 00256 */ 00257 # include <cppad/local/std_set.hpp> 00258 # include <cppad/local/color_general.hpp> 00259 00260 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00261 /*! 00262 \file sparse_jacobian.hpp 00263 Sparse Jacobian driver routine and helper functions. 00264 */ 00265 // =========================================================================== 00266 /*! 00267 class used by SparseJacobian to hold information so it does not need to be 00268 recomputed. 00269 */ 00270 class sparse_jacobian_work { 00271 public: 00272 /// Coloring method: color_general_cppad or color_general_colpack 00273 /// (this field is set by user) 00274 std::string color_method; 00275 /// indices that sort the user row and col arrays by color 00276 CppAD::vector<size_t> order; 00277 /// results of the coloring algorithm 00278 CppAD::vector<size_t> color; 00279 00280 /// constructor 00281 sparse_jacobian_work(void) : color_method("cppad") 00282 { } 00283 /// reset coloring method to its default and 00284 /// inform CppAD that color and order need to be recomputed 00285 void clear(void) 00286 { color_method = "cppad"; 00287 order.clear(); 00288 color.clear(); 00289 } 00290 }; 00291 // =========================================================================== 00292 /*! 00293 Private helper function forward mode cases 00294 00295 \tparam Base 00296 is the base type for the recording that is stored in this 00297 <code>ADFun<Base></code> object. 00298 00299 \tparam VectorBase 00300 is a simple vector class with elements of type \a Base. 00301 00302 \tparam VectorSet 00303 is either \c sparse_pack, \c sparse_set or \c sparse_list. 00304 00305 \tparam VectorSize 00306 is a simple vector class with elements of type \c size_t. 00307 00308 \param x [in] 00309 is a vector specifing the point at which to compute the Jacobian. 00310 00311 \param p_transpose [in] 00312 If <code>work.color.size() != 0</code>, 00313 then \c p_transpose is not used. 00314 Otherwise, it is a 00315 sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object. 00316 Note that we do not change the values in \c p_transpose, 00317 but is not \c const because we use its iterator facility. 00318 00319 \param row [in] 00320 is the vector of row indices for the returned Jacobian values. 00321 00322 \param col [in] 00323 is the vector of columns indices for the returned Jacobian values. 00324 It must have the same size as \c row. 00325 00326 \param jac [out] 00327 is the vector of Jacobian values. We use \c K to denote the size of \c jac. 00328 The return value <code>jac[k]</code> is the partial of the 00329 <code>row[k]</code> range component of the function with respect 00330 the the <code>col[k]</code> domain component of its argument. 00331 00332 \param work 00333 <code>work.color_method</code> is an input. The rest of 00334 this structure contains information that is computed by \c SparseJacobainFor. 00335 If the sparsity pattern, \c row vector, or \c col vectors 00336 are not the same between calls to \c SparseJacobianFor, 00337 \c work.clear() must be called to reinitialize \c work. 00338 00339 \return 00340 Is the number of first order forward sweeps used to compute the 00341 requested Jacobian values. The total work, not counting the zero order 00342 forward sweep, or the time to combine computations, is proportional to this 00343 return value. 00344 */ 00345 template<class Base> 00346 template <class VectorBase, class VectorSet, class VectorSize> 00347 size_t ADFun<Base>::SparseJacobianFor( 00348 const VectorBase& x , 00349 VectorSet& p_transpose , 00350 const VectorSize& row , 00351 const VectorSize& col , 00352 VectorBase& jac , 00353 sparse_jacobian_work& work ) 00354 { 00355 size_t j, k, ell; 00356 00357 CppAD::vector<size_t>& order(work.order); 00358 CppAD::vector<size_t>& color(work.color); 00359 00360 size_t m = Range(); 00361 size_t n = Domain(); 00362 00363 // some values 00364 const Base zero(0); 00365 const Base one(1); 00366 00367 // check VectorBase is Simple Vector class with Base type elements 00368 CheckSimpleVector<Base, VectorBase>(); 00369 00370 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00371 CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); 00372 00373 // number of components of Jacobian that are required 00374 size_t K = size_t(jac.size()); 00375 CPPAD_ASSERT_UNKNOWN( row.size() == K ); 00376 CPPAD_ASSERT_UNKNOWN( col.size() == K ); 00377 00378 // Point at which we are evaluating the Jacobian 00379 Forward(0, x); 00380 00381 // check for case where nothing (except Forward above) to do 00382 if( K == 0 ) 00383 return 0; 00384 00385 if( color.size() == 0 ) 00386 { 00387 CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() == n ); 00388 CPPAD_ASSERT_UNKNOWN( p_transpose.end() == m ); 00389 00390 // execute coloring algorithm 00391 color.resize(n); 00392 if( work.color_method == "cppad" ) 00393 color_general_cppad(p_transpose, col, row, color); 00394 else if( work.color_method == "colpack" ) 00395 { 00396 # if CPPAD_HAS_COLPACK 00397 color_general_colpack(p_transpose, col, row, color); 00398 # else 00399 CPPAD_ASSERT_KNOWN( 00400 false, 00401 "SparseJacobianForward: work.color_method = colpack " 00402 "and colpack_prefix missing from cmake command line." 00403 ); 00404 # endif 00405 } 00406 else CPPAD_ASSERT_KNOWN( 00407 false, 00408 "SparseJacobianForward: work.color_method is not valid." 00409 ); 00410 00411 // put sorting indices in color order 00412 VectorSize key(K); 00413 order.resize(K); 00414 for(k = 0; k < K; k++) 00415 key[k] = color[ col[k] ]; 00416 index_sort(key, order); 00417 } 00418 size_t n_color = 1; 00419 for(j = 0; j < n; j++) if( color[j] < n ) 00420 n_color = std::max(n_color, color[j] + 1); 00421 00422 // direction vector for calls to forward 00423 VectorBase dx(n); 00424 00425 // location for return values from forward 00426 VectorBase dy(m); 00427 00428 // initialize the return value 00429 for(k = 0; k < K; k++) 00430 jac[k] = zero; 00431 00432 // loop over colors 00433 k = 0; 00434 for(ell = 0; ell < n_color; ell++) 00435 { CPPAD_ASSERT_UNKNOWN( color[ col[ order[k] ] ] == ell ); 00436 00437 // combine all columns with this color 00438 for(j = 0; j < n; j++) 00439 { dx[j] = zero; 00440 if( color[j] == ell ) 00441 dx[j] = one; 00442 } 00443 // call forward mode for all these columns at once 00444 dy = Forward(1, dx); 00445 00446 // set the corresponding components of the result 00447 while( k < K && color[ col[order[k]] ] == ell ) 00448 { jac[ order[k] ] = dy[row[order[k]]]; 00449 k++; 00450 } 00451 } 00452 return n_color; 00453 } 00454 /*! 00455 Private helper function for reverse mode cases. 00456 00457 \tparam Base 00458 is the base type for the recording that is stored in this 00459 <code>ADFun<Base></code> object. 00460 00461 \tparam VectorBase 00462 is a simple vector class with elements of type \a Base. 00463 00464 \tparam VectorSet 00465 is either \c sparse_pack, \c sparse_set or \c sparse_list. 00466 00467 \tparam VectorSize 00468 is a simple vector class with elements of type \c size_t. 00469 00470 \param x [in] 00471 is a vector specifing the point at which to compute the Jacobian. 00472 00473 \param p [in] 00474 If <code>work.color.size() != 0</code>, then \c p is not used. 00475 Otherwise, it is a 00476 sparsity pattern for the Jacobian of this ADFun<Base> object. 00477 Note that we do not change the values in \c p, 00478 but is not \c const because we use its iterator facility. 00479 00480 \param row [in] 00481 is the vector of row indices for the returned Jacobian values. 00482 00483 \param col [in] 00484 is the vector of columns indices for the returned Jacobian values. 00485 It must have the same size as \c row. 00486 00487 \param jac [out] 00488 is the vector of Jacobian values. 00489 It must have the same size as \c row. 00490 The return value <code>jac[k]</code> is the partial of the 00491 <code>row[k]</code> range component of the function with respect 00492 the the <code>col[k]</code> domain component of its argument. 00493 00494 \param work 00495 <code>work.color_method</code> is an input. The rest of 00496 This structure contains information that is computed by \c SparseJacobainRev. 00497 If the sparsity pattern, \c row vector, or \c col vectors 00498 are not the same between calls to \c SparseJacobianRev, 00499 \c work.clear() must be called to reinitialize \c work. 00500 00501 \return 00502 Is the number of first order reverse sweeps used to compute the 00503 reverse Jacobian values. The total work, not counting the zero order 00504 forward sweep, or the time to combine computations, is proportional to this 00505 return value. 00506 */ 00507 template<class Base> 00508 template <class VectorBase, class VectorSet, class VectorSize> 00509 size_t ADFun<Base>::SparseJacobianRev( 00510 const VectorBase& x , 00511 VectorSet& p , 00512 const VectorSize& row , 00513 const VectorSize& col , 00514 VectorBase& jac , 00515 sparse_jacobian_work& work ) 00516 { 00517 size_t i, k, ell; 00518 00519 CppAD::vector<size_t>& order(work.order); 00520 CppAD::vector<size_t>& color(work.color); 00521 00522 size_t m = Range(); 00523 size_t n = Domain(); 00524 00525 // some values 00526 const Base zero(0); 00527 const Base one(1); 00528 00529 // check VectorBase is Simple Vector class with Base type elements 00530 CheckSimpleVector<Base, VectorBase>(); 00531 00532 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00533 CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 ); 00534 00535 // number of components of Jacobian that are required 00536 size_t K = size_t(jac.size()); 00537 CPPAD_ASSERT_UNKNOWN( row.size() == K ); 00538 CPPAD_ASSERT_UNKNOWN( col.size() == K ); 00539 00540 // Point at which we are evaluating the Jacobian 00541 Forward(0, x); 00542 00543 // check for case where nothing (except Forward above) to do 00544 if( K == 0 ) 00545 return 0; 00546 00547 if( color.size() == 0 ) 00548 { 00549 CPPAD_ASSERT_UNKNOWN( p.n_set() == m ); 00550 CPPAD_ASSERT_UNKNOWN( p.end() == n ); 00551 00552 // execute the coloring algorithm 00553 color.resize(m); 00554 if( work.color_method == "cppad" ) 00555 color_general_cppad(p, row, col, color); 00556 else if( work.color_method == "colpack" ) 00557 { 00558 # if CPPAD_HAS_COLPACK 00559 color_general_colpack(p, row, col, color); 00560 # else 00561 CPPAD_ASSERT_KNOWN( 00562 false, 00563 "SparseJacobianReverse: work.color_method = colpack " 00564 "and colpack_prefix missing from cmake command line." 00565 ); 00566 # endif 00567 } 00568 else CPPAD_ASSERT_KNOWN( 00569 false, 00570 "SparseJacobianReverse: work.color_method is not valid." 00571 ); 00572 00573 // put sorting indices in color order 00574 VectorSize key(K); 00575 order.resize(K); 00576 for(k = 0; k < K; k++) 00577 key[k] = color[ row[k] ]; 00578 index_sort(key, order); 00579 } 00580 size_t n_color = 1; 00581 for(i = 0; i < m; i++) if( color[i] < m ) 00582 n_color = std::max(n_color, color[i] + 1); 00583 00584 // weighting vector for calls to reverse 00585 VectorBase w(m); 00586 00587 // location for return values from Reverse 00588 VectorBase dw(n); 00589 00590 // initialize the return value 00591 for(k = 0; k < K; k++) 00592 jac[k] = zero; 00593 00594 // loop over colors 00595 k = 0; 00596 for(ell = 0; ell < n_color; ell++) 00597 { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell ); 00598 00599 // combine all the rows with this color 00600 for(i = 0; i < m; i++) 00601 { w[i] = zero; 00602 if( color[i] == ell ) 00603 w[i] = one; 00604 } 00605 // call reverse mode for all these rows at once 00606 dw = Reverse(1, w); 00607 00608 // set the corresponding components of the result 00609 while( k < K && color[ row[order[k]] ] == ell ) 00610 { jac[ order[k] ] = dw[col[order[k]]]; 00611 k++; 00612 } 00613 } 00614 return n_color; 00615 } 00616 // ========================================================================== 00617 // Public Member functions 00618 // ========================================================================== 00619 /*! 00620 Compute user specified subset of a sparse Jacobian using forward mode. 00621 00622 The C++ source code corresponding to this operation is 00623 \verbatim 00624 SparceJacobianForward(x, p, row, col, jac, work) 00625 \endverbatim 00626 00627 \tparam Base 00628 is the base type for the recording that is stored in this 00629 <code>ADFun<Base></code> object. 00630 00631 \tparam VectorBase 00632 is a simple vector class with elements of type \a Base. 00633 00634 \tparam VectorSet 00635 is a simple vector class with elements of type 00636 \c bool or \c std::set<size_t>. 00637 00638 \tparam VectorSize 00639 is a simple vector class with elements of type \c size_t. 00640 00641 \param x [in] 00642 is a vector specifing the point at which to compute the Jacobian. 00643 00644 \param p [in] 00645 is the sparsity pattern for the Jacobian that we are calculating. 00646 00647 \param row [in] 00648 is the vector of row indices for the returned Jacobian values. 00649 00650 \param col [in] 00651 is the vector of columns indices for the returned Jacobian values. 00652 It must have the same size as \c row. 00653 00654 \param jac [out] 00655 is the vector of Jacobian values. 00656 It must have the same size as \c row. 00657 The return value <code>jac[k]</code> is the partial of the 00658 <code>row[k]</code> range component of the function with respect 00659 the the <code>col[k]</code> domain component of its argument. 00660 00661 \param work [in,out] 00662 this structure contains information that depends on the function object, 00663 sparsity pattern, \c row vector, and \c col vector. 00664 If they are not the same between calls to \c SparseJacobianForward, 00665 \c work.clear() must be called to reinitialize them. 00666 00667 \return 00668 Is the number of first order forward sweeps used to compute the 00669 requested Jacobian values. The total work, not counting the zero order 00670 forward sweep, or the time to combine computations, is proportional to this 00671 return value. 00672 */ 00673 template<class Base> 00674 template <class VectorBase, class VectorSet, class VectorSize> 00675 size_t ADFun<Base>::SparseJacobianForward( 00676 const VectorBase& x , 00677 const VectorSet& p , 00678 const VectorSize& row , 00679 const VectorSize& col , 00680 VectorBase& jac , 00681 sparse_jacobian_work& work ) 00682 { 00683 size_t n = Domain(); 00684 size_t m = Range(); 00685 # ifndef NDEBUG 00686 size_t k, K = jac.size(); 00687 CPPAD_ASSERT_KNOWN( 00688 size_t(x.size()) == n , 00689 "SparseJacobianForward: size of x not equal domain dimension for f." 00690 ); 00691 CPPAD_ASSERT_KNOWN( 00692 size_t(row.size()) == K && size_t(col.size()) == K , 00693 "SparseJacobianForward: either r or c does not have " 00694 "the same size as jac." 00695 ); 00696 CPPAD_ASSERT_KNOWN( 00697 work.color.size() == 0 || work.color.size() == n, 00698 "SparseJacobianForward: invalid value in work." 00699 ); 00700 for(k = 0; k < K; k++) 00701 { CPPAD_ASSERT_KNOWN( 00702 row[k] < m, 00703 "SparseJacobianForward: invalid value in r." 00704 ); 00705 CPPAD_ASSERT_KNOWN( 00706 col[k] < n, 00707 "SparseJacobianForward: invalid value in c." 00708 ); 00709 } 00710 if( work.color.size() != 0 ) 00711 for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN( 00712 work.color[j] <= n, 00713 "SparseJacobianForward: invalid value in work." 00714 ); 00715 # endif 00716 00717 typedef typename VectorSet::value_type Set_type; 00718 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00719 Pattern_type s_transpose; 00720 if( work.color.size() == 0 ) 00721 { bool transpose = true; 00722 sparsity_user2internal(s_transpose, p, m, n, transpose); 00723 } 00724 size_t n_sweep = SparseJacobianFor(x, s_transpose, row, col, jac, work); 00725 return n_sweep; 00726 } 00727 /*! 00728 Compute user specified subset of a sparse Jacobian using forward mode. 00729 00730 The C++ source code corresponding to this operation is 00731 \verbatim 00732 SparceJacobianReverse(x, p, row, col, jac, work) 00733 \endverbatim 00734 00735 \tparam Base 00736 is the base type for the recording that is stored in this 00737 <code>ADFun<Base></code> object. 00738 00739 \tparam VectorBase 00740 is a simple vector class with elements of type \a Base. 00741 00742 \tparam VectorSet 00743 is a simple vector class with elements of type 00744 \c bool or \c std::set<size_t>. 00745 00746 \tparam VectorSize 00747 is a simple vector class with elements of type \c size_t. 00748 00749 \param x [in] 00750 is a vector specifing the point at which to compute the Jacobian. 00751 00752 \param p [in] 00753 is the sparsity pattern for the Jacobian that we are calculating. 00754 00755 \param row [in] 00756 is the vector of row indices for the returned Jacobian values. 00757 00758 \param col [in] 00759 is the vector of columns indices for the returned Jacobian values. 00760 It must have the same size as \c row. 00761 00762 \param jac [out] 00763 is the vector of Jacobian values. 00764 It must have the same size as \c row. 00765 The return value <code>jac[k]</code> is the partial of the 00766 <code>row[k]</code> range component of the function with respect 00767 the the <code>col[k]</code> domain component of its argument. 00768 00769 \param work [in,out] 00770 this structure contains information that depends on the function object, 00771 sparsity pattern, \c row vector, and \c col vector. 00772 If they are not the same between calls to \c SparseJacobianReverse, 00773 \c work.clear() must be called to reinitialize them. 00774 00775 \return 00776 Is the number of first order reverse sweeps used to compute the 00777 reverse Jacobian values. The total work, not counting the zero order 00778 forward sweep, or the time to combine computations, is proportional to this 00779 return value. 00780 */ 00781 template<class Base> 00782 template <class VectorBase, class VectorSet, class VectorSize> 00783 size_t ADFun<Base>::SparseJacobianReverse( 00784 const VectorBase& x , 00785 const VectorSet& p , 00786 const VectorSize& row , 00787 const VectorSize& col , 00788 VectorBase& jac , 00789 sparse_jacobian_work& work ) 00790 { 00791 size_t m = Range(); 00792 size_t n = Domain(); 00793 # ifndef NDEBUG 00794 size_t k, K = jac.size(); 00795 CPPAD_ASSERT_KNOWN( 00796 size_t(x.size()) == n , 00797 "SparseJacobianReverse: size of x not equal domain dimension for f." 00798 ); 00799 CPPAD_ASSERT_KNOWN( 00800 size_t(row.size()) == K && size_t(col.size()) == K , 00801 "SparseJacobianReverse: either r or c does not have " 00802 "the same size as jac." 00803 ); 00804 CPPAD_ASSERT_KNOWN( 00805 work.color.size() == 0 || work.color.size() == m, 00806 "SparseJacobianReverse: invalid value in work." 00807 ); 00808 for(k = 0; k < K; k++) 00809 { CPPAD_ASSERT_KNOWN( 00810 row[k] < m, 00811 "SparseJacobianReverse: invalid value in r." 00812 ); 00813 CPPAD_ASSERT_KNOWN( 00814 col[k] < n, 00815 "SparseJacobianReverse: invalid value in c." 00816 ); 00817 } 00818 if( work.color.size() != 0 ) 00819 for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 00820 work.color[i] <= m, 00821 "SparseJacobianReverse: invalid value in work." 00822 ); 00823 # endif 00824 00825 typedef typename VectorSet::value_type Set_type; 00826 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00827 Pattern_type s; 00828 if( work.color.size() == 0 ) 00829 { bool transpose = false; 00830 sparsity_user2internal(s, p, m, n, transpose); 00831 } 00832 size_t n_sweep = SparseJacobianRev(x, s, row, col, jac, work); 00833 return n_sweep; 00834 } 00835 /*! 00836 Compute a sparse Jacobian. 00837 00838 The C++ source code corresponding to this operation is 00839 \verbatim 00840 jac = SparseJacobian(x, p) 00841 \endverbatim 00842 00843 \tparam Base 00844 is the base type for the recording that is stored in this 00845 <code>ADFun<Base></code> object. 00846 00847 \tparam VectorBase 00848 is a simple vector class with elements of type \a Base. 00849 00850 \tparam VectorSet 00851 is a simple vector class with elements of type 00852 \c bool or \c std::set<size_t>. 00853 00854 \param x [in] 00855 is a vector specifing the point at which to compute the Jacobian. 00856 00857 \param p [in] 00858 is the sparsity pattern for the Jacobian that we are calculating. 00859 00860 \return 00861 Will be a vector if size \c m * n containing the Jacobian at the 00862 specified point (in row major order). 00863 */ 00864 template <class Base> 00865 template <class VectorBase, class VectorSet> 00866 VectorBase ADFun<Base>::SparseJacobian( 00867 const VectorBase& x, const VectorSet& p 00868 ) 00869 { size_t i, j, k; 00870 00871 size_t m = Range(); 00872 size_t n = Domain(); 00873 VectorBase jac(m * n); 00874 00875 CPPAD_ASSERT_KNOWN( 00876 size_t(x.size()) == n, 00877 "SparseJacobian: size of x not equal domain size for f." 00878 ); 00879 CheckSimpleVector<Base, VectorBase>(); 00880 00881 typedef typename VectorSet::value_type Set_type; 00882 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00883 00884 // initialize the return value as zero 00885 Base zero(0); 00886 for(i = 0; i < m; i++) 00887 for(j = 0; j < n; j++) 00888 jac[i * n + j] = zero; 00889 00890 sparse_jacobian_work work; 00891 CppAD::vector<size_t> row; 00892 CppAD::vector<size_t> col; 00893 if( n <= m ) 00894 { 00895 // need an internal copy of sparsity pattern 00896 Pattern_type s_transpose; 00897 bool transpose = true; 00898 sparsity_user2internal(s_transpose, p, m, n, transpose); 00899 00900 k = 0; 00901 for(j = 0; j < n; j++) 00902 { s_transpose.begin(j); 00903 i = s_transpose.next_element(); 00904 while( i != s_transpose.end() ) 00905 { row.push_back(i); 00906 col.push_back(j); 00907 k++; 00908 i = s_transpose.next_element(); 00909 } 00910 } 00911 size_t K = k; 00912 VectorBase J(K); 00913 00914 // now we have folded this into the following case 00915 SparseJacobianFor(x, s_transpose, row, col, J, work); 00916 00917 // now set the non-zero return values 00918 for(k = 0; k < K; k++) 00919 jac[ row[k] * n + col[k] ] = J[k]; 00920 } 00921 else 00922 { 00923 // need an internal copy of sparsity pattern 00924 Pattern_type s; 00925 bool transpose = false; 00926 sparsity_user2internal(s, p, m, n, transpose); 00927 00928 k = 0; 00929 for(i = 0; i < m; i++) 00930 { s.begin(i); 00931 j = s.next_element(); 00932 while( j != s.end() ) 00933 { row.push_back(i); 00934 col.push_back(j); 00935 k++; 00936 j = s.next_element(); 00937 } 00938 } 00939 size_t K = k; 00940 VectorBase J(K); 00941 00942 // now we have folded this into the following case 00943 SparseJacobianRev(x, s, row, col, J, work); 00944 00945 // now set the non-zero return values 00946 for(k = 0; k < K; k++) 00947 jac[ row[k] * n + col[k] ] = J[k]; 00948 } 00949 00950 return jac; 00951 } 00952 00953 /*! 00954 Compute a sparse Jacobian. 00955 00956 The C++ source code corresponding to this operation is 00957 \verbatim 00958 jac = SparseJacobian(x) 00959 \endverbatim 00960 00961 \tparam Base 00962 is the base type for the recording that is stored in this 00963 <code>ADFun<Base></code> object. 00964 00965 \tparam VectorBase 00966 is a simple vector class with elements of the \a Base. 00967 00968 \param x [in] 00969 is a vector specifing the point at which to compute the Jacobian. 00970 00971 \return 00972 Will be a vector of size \c m * n containing the Jacobian at the 00973 specified point (in row major order). 00974 */ 00975 template <class Base> 00976 template <class VectorBase> 00977 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x ) 00978 { typedef CppAD::vectorBool VectorBool; 00979 00980 size_t m = Range(); 00981 size_t n = Domain(); 00982 00983 // sparsity pattern for Jacobian 00984 VectorBool p(m * n); 00985 00986 if( n <= m ) 00987 { size_t j, k; 00988 00989 // use forward mode 00990 VectorBool r(n * n); 00991 for(j = 0; j < n; j++) 00992 { for(k = 0; k < n; k++) 00993 r[j * n + k] = false; 00994 r[j * n + j] = true; 00995 } 00996 p = ForSparseJac(n, r); 00997 } 00998 else 00999 { size_t i, k; 01000 01001 // use reverse mode 01002 VectorBool s(m * m); 01003 for(i = 0; i < m; i++) 01004 { for(k = 0; k < m; k++) 01005 s[i * m + k] = false; 01006 s[i * m + i] = true; 01007 } 01008 p = RevSparseJac(m, s); 01009 } 01010 return SparseJacobian(x, p); 01011 } 01012 01013 } // END_CPPAD_NAMESPACE 01014 # endif