CppAD: A C++ Algorithmic Differentiation Package
20130102
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_JACOBIAN_INCLUDED 00003 # define CPPAD_SPARSE_JACOBIAN_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 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 recomputed 00019 valarray 00020 std 00021 CppAD 00022 Bool 00023 jac 00024 Jacobian 00025 const 00026 Taylor 00027 $$ 00028 00029 $section Sparse Jacobian: Easy Driver$$ 00030 $index SparseJacobian$$ 00031 $index jacobian, sparse$$ 00032 00033 $head Syntax$$ 00034 $icode%jac% = %f%.SparseJacobian(%x%) 00035 %jac% = %f%.SparseJacobian(%x%, %p%) 00036 %n_sweep% = %f%.SparseJacobianForward(%x%, %p%, %row%, %col%, %jac%, %work%) 00037 %n_sweep% = %f%.SparseJacobianReverse(%x%, %p%, %row%, %col%, %jac%, %work%) 00038 %$$ 00039 00040 $head Purpose$$ 00041 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size, 00042 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$. 00043 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the 00044 $cref/AD function/glossary/AD Function/$$ 00045 corresponding to $icode f$$. 00046 The syntax above sets $icode jac$$ to the Jacobian 00047 $latex \[ 00048 jac = F^{(1)} (x) 00049 \] $$ 00050 This routine takes advantage of the sparsity of the Jacobian 00051 in order to reduce the amount of computation necessary. 00052 If $icode row$$ and $icode col$$ are present, it also takes 00053 advantage of the reduced set of elements of the Jacobian that 00054 need to be computed. 00055 One can use speed tests (e.g. $cref speed_test$$) 00056 to verify that results are computed faster 00057 than when using the routine $cref Jacobian$$. 00058 00059 $head f$$ 00060 The object $icode f$$ has prototype 00061 $codei% 00062 ADFun<%Base%> %f% 00063 %$$ 00064 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$ 00065 (see $cref/Uses Forward/sparse_jacobian/Uses Forward/$$ below). 00066 00067 $head x$$ 00068 The argument $icode x$$ has prototype 00069 $codei% 00070 const %VectorBase%& %x% 00071 %$$ 00072 (see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below) 00073 and its size 00074 must be equal to $icode n$$, the dimension of the 00075 $cref/domain/seq_property/Domain/$$ space for $icode f$$. 00076 It specifies 00077 that point at which to evaluate the Jacobian. 00078 00079 $head p$$ 00080 The argument $icode p$$ is optional and has prototype 00081 $codei% 00082 const %VectorSet%& %p% 00083 %$$ 00084 (see $cref/VectorSet/sparse_jacobian/VectorSet/$$ below). 00085 If it has elements of type $code bool$$, 00086 its size is $latex m * n$$. 00087 If it has elements of type $code std::set<size_t>$$, 00088 its size is $latex m$$ and all its set elements are between 00089 zero and $latex n - 1$$. 00090 It specifies a 00091 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00092 for the Jacobian $latex F^{(1)} (x)$$. 00093 $pre 00094 00095 $$ 00096 If this sparsity pattern does not change between calls to 00097 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 00098 (using $cref ForSparseJac$$ or $cref RevSparseJac$$) 00099 and then pass $icode p$$ to $codei SparseJacobian$$. 00100 In addition, 00101 if you specify $icode p$$, CppAD will use the same 00102 type of sparsity representation 00103 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00104 for its internal calculations. 00105 Otherwise, the representation 00106 for the internal calculations is unspecified. 00107 00108 $head row, col$$ 00109 The arguments $icode row$$ and $icode col$$ are optional and have prototype 00110 $codei% 00111 const %VectorSize%& %row% 00112 const %VectorSize%& %col% 00113 %$$ 00114 (see $cref/VectorSize/sparse_jacobian/VectorSize/$$ below). 00115 They specify which rows and columns of $latex F^{(1)} (x)$$ are 00116 returned and in what order. 00117 We use $latex K$$ to denote the value $icode%jac%.size()%$$ 00118 which must also equal the size of $icode row$$ and $icode col$$. 00119 Furthermore, 00120 for $latex k = 0 , \ldots , K-1$$, it must hold that 00121 $latex row[k] < m$$ and $latex col[k] < n$$. 00122 In addition, 00123 all of the $latex (row[k], col[k])$$ pairs must correspond to a true value 00124 in the sparsity pattern $icode p$$. 00125 00126 $head jac$$ 00127 The result $icode jac$$ has prototype 00128 $codei% 00129 %VectorBase%& %jac% 00130 %$$ 00131 In the case where the arguments $icode row$$ and $icode col$$ are not present, 00132 the size of $icode jac$$ is $latex m * n$$ and 00133 for $latex i = 0 , \ldots , m-1$$, 00134 $latex j = 0 , \ldots , n-1$$, 00135 $latex \[ 00136 jac [ i * n + j ] = \D{ F_i }{ x_j } (x) 00137 \] $$ 00138 $pre 00139 00140 $$ 00141 In the case where the arguments $icode row$$ and $icode col$$ are present, 00142 we use $latex K$$ to denote the size of $icode jac$$. 00143 The input value of its elements does not matter. 00144 Upon return, for $latex k = 0 , \ldots , K - 1$$, 00145 $latex \[ 00146 jac [ k ] = \D{ F_i }{ x_j } (x) 00147 \; , \; 00148 \; {\rm where} \; 00149 i = row[k] 00150 \; {\rm and } \; 00151 j = col[k] 00152 \] $$ 00153 00154 $head work$$ 00155 If this argument is present, it has prototype 00156 $icode% 00157 sparse_jacobian_work& %work% 00158 %$$ 00159 This object can only be used with the routines 00160 $code SparseJacobianForward$$ and $code SparseJacobianReverse$$. 00161 During its the first use, information is stored in $icode work$$. 00162 This is used to reduce the work done by future calls to the same mode 00163 (forward or reverse), 00164 the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$. 00165 If a future call is for a different mode, 00166 or any of these values have changed, 00167 you must first call $icode%work%.clear()%$$ 00168 to inform CppAD that this information needs to be recomputed. 00169 00170 $head n_sweep$$ 00171 The return value $icode n_sweep$$ has prototype 00172 $codei% 00173 size_t %n_sweep% 00174 %$$ 00175 If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used, 00176 $icode n_sweep$$ is the number of first order forward (reverse) sweeps 00177 used to compute the requested Jacobian values. 00178 This is proportional to the total work that $code SparseJacobian$$ does, 00179 not counting the zero order forward sweep, 00180 or the work to combine multiple columns (rows) into a single sweep. 00181 00182 $head VectorBase$$ 00183 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with 00184 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00185 $icode Base$$. 00186 The routine $cref CheckSimpleVector$$ will generate an error message 00187 if this is not the case. 00188 00189 $head VectorSet$$ 00190 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with 00191 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00192 $code bool$$ or $code std::set<size_t>$$; 00193 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00194 of the difference. 00195 The routine $cref CheckSimpleVector$$ will generate an error message 00196 if this is not the case. 00197 00198 $subhead Restrictions$$ 00199 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00200 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00201 corresponding set. 00202 According to section 26.3.2.3 of the 1998 C++ standard, 00203 $code std::valarray< std::set<size_t> >$$ does not satisfy 00204 this condition. 00205 00206 $head VectorSize$$ 00207 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with 00208 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00209 $code size_t$$. 00210 The routine $cref CheckSimpleVector$$ will generate an error message 00211 if this is not the case. 00212 00213 $head Uses Forward$$ 00214 After each call to $cref Forward$$, 00215 the object $icode f$$ contains the corresponding 00216 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00217 After a call to any of the sparse Jacobian routines, 00218 the zero order Taylor coefficients correspond to 00219 $icode%f%.Forward(0, %x%)%$$ 00220 and the other coefficients are unspecified. 00221 00222 After $code SparseJacobian$$, 00223 the previous calls to $cref Forward$$ are undefined. 00224 00225 $head Example$$ 00226 $children% 00227 example/sparse_jacobian.cpp 00228 %$$ 00229 The routine 00230 $cref sparse_jacobian.cpp$$ 00231 is examples and tests of $code sparse_jacobian$$. 00232 It return $code true$$, if it succeeds and $code false$$ otherwise. 00233 00234 $end 00235 ============================================================================== 00236 */ 00237 # include <cppad/local/std_set.hpp> 00238 00239 CPPAD_BEGIN_NAMESPACE 00240 /*! 00241 \defgroup sparse_jacobian_hpp sparse_jacobian.hpp 00242 \{ 00243 \file sparse_jacobian.hpp 00244 Sparse Jacobian driver routine and helper functions. 00245 */ 00246 // =========================================================================== 00247 /*! 00248 class used by SparseJacobian to hold information so it does not need to be 00249 recomputed. 00250 */ 00251 class sparse_jacobian_work { 00252 public: 00253 /// version of user row array with the extra value m at end 00254 CppAD::vector<size_t> user_row; 00255 /// version of user col array with the extra value n at end 00256 CppAD::vector<size_t> user_col; 00257 /// indices that sort the user arrays by row 00258 /// with the extra value K at the end 00259 CppAD::vector<size_t> sort_row; 00260 /// indices that sort the user arrays by column 00261 /// with the extra value K at the end 00262 CppAD::vector<size_t> sort_col; 00263 /// results of the coloring algorithm 00264 CppAD::vector<size_t> color; 00265 /// inform CppAD that this information needs to be recomputed 00266 void clear(void) 00267 { user_row.clear(); 00268 user_col.clear(); 00269 sort_row.clear(); 00270 sort_col.clear(); 00271 color.clear(); 00272 } 00273 }; 00274 // =========================================================================== 00275 /*! 00276 Private helper function forward mode cases 00277 00278 \tparam Base 00279 See \c SparseJacobianForward(x, p, row, col, jac, work). 00280 00281 \tparam VectorBase 00282 See \c SparseJacobianForward(x, p, row, col, jac, work). 00283 00284 \tparam VectorSet 00285 is either \c sparse_pack, \c sparse_set or \c sparse_list. 00286 00287 \param x 00288 See \c SparseJacobianForward(x, p, row, col, jac, work). 00289 00290 \param p_transpose 00291 If <code>work.color.size() != 0</code>, 00292 then \c p_transpose is not used. 00293 Otherwise, it is a 00294 sparsity pattern for the transpose of the Jacobian of this ADFun<Base> object. 00295 Note that we do not change the values in \c p_transpose, 00296 but is not \c const because we use its iterator facility. 00297 00298 00299 \param jac 00300 See \c SparseJacobianForward(x, p, row, col, jac, work). 00301 00302 \param work 00303 See \c SparseJacobianForward(x, p, row, col, jac, work). 00304 In addition, \c work.color contains information that is computed by 00305 \c SparseJacobainFor. 00306 If the sparsity pattern, \c row vector, and \c col vectors 00307 are not the same between calls, \c work.clear() must be called 00308 to reinitialize \c work.color. 00309 00310 \return 00311 See \c SparseJacobianForward(x, p, row, col, jac, work). 00312 */ 00313 template<class Base> 00314 template <class VectorBase, class VectorSet> 00315 size_t ADFun<Base>::SparseJacobianFor( 00316 const VectorBase& x , 00317 VectorSet& p_transpose , 00318 VectorBase& jac , 00319 sparse_jacobian_work& work ) 00320 { 00321 using CppAD::vectorBool; 00322 size_t i, j, k, ell; 00323 00324 CppAD::vector<size_t>& row(work.user_row); 00325 CppAD::vector<size_t>& col(work.user_col); 00326 CppAD::vector<size_t>& sort_col(work.sort_col); 00327 CppAD::vector<size_t>& color(work.color); 00328 00329 size_t m = Range(); 00330 size_t n = Domain(); 00331 00332 // some values 00333 const Base zero(0); 00334 const Base one(1); 00335 00336 // check VectorBase is Simple Vector class with Base type elements 00337 CheckSimpleVector<Base, VectorBase>(); 00338 00339 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00340 CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); 00341 00342 // number of components of Jacobian that are required 00343 size_t K = size_t(jac.size()); 00344 CPPAD_ASSERT_UNKNOWN( row.size() == K+1 ); 00345 CPPAD_ASSERT_UNKNOWN( col.size() == K+1 ); 00346 CPPAD_ASSERT_UNKNOWN( row[K] == m ); 00347 CPPAD_ASSERT_UNKNOWN( col[K] == n ); 00348 00349 // Point at which we are evaluating the Jacobian 00350 Forward(0, x); 00351 00352 if( color.size() == 0 ) 00353 { 00354 CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() == n ); 00355 CPPAD_ASSERT_UNKNOWN( p_transpose.end() == m ); 00356 00357 // rows and columns that are in the returned jacobian 00358 VectorSet r_used, c_used; 00359 r_used.resize(n, m); 00360 c_used.resize(m, n); 00361 k = 0; 00362 while( k < K ) 00363 { CPPAD_ASSERT_UNKNOWN( 00364 row[sort_col[k]] < m && col[sort_col[k]] < n 00365 ); 00366 CPPAD_ASSERT_UNKNOWN( 00367 k == 0 || col[sort_col[k-1]] <= col[sort_col[k]] 00368 ); 00369 CPPAD_ASSERT_KNOWN( 00370 p_transpose.is_element(col[sort_col[k]], row[sort_col[k]]) , 00371 "SparseJacobianForward: " 00372 "an (row, col) pair is not in sparsity pattern." 00373 ); 00374 r_used.add_element(col[sort_col[k]], row[sort_col[k]]); 00375 c_used.add_element(row[sort_col[k]], col[sort_col[k]]); 00376 k++; 00377 } 00378 00379 // given a row index, which columns are non-zero and not used 00380 VectorSet not_used; 00381 not_used.resize(m, n); 00382 for(j = 0; j < n; j++) 00383 { p_transpose.begin(j); 00384 i = p_transpose.next_element(); 00385 while( i != p_transpose.end() ) 00386 { if( ! c_used.is_element(i, j) ) 00387 not_used.add_element(i, j); 00388 i = p_transpose.next_element(); 00389 } 00390 } 00391 00392 // initial coloring 00393 color.resize(n); 00394 for(j = 0; j < n; j++) 00395 color[j] = j; 00396 00397 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00398 // Graph Coloring in Optimization Revisited by 00399 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00400 vectorBool forbidden(n); 00401 for(j = 1; j < n; j++) 00402 { 00403 // initial all colors as ok for this column 00404 // (value of forbidden for ell > j does not matter) 00405 for(ell = 0; ell <= j; ell++) 00406 forbidden[ell] = false; 00407 00408 // for each row that is non-zero for this column 00409 p_transpose.begin(j); 00410 i = p_transpose.next_element(); 00411 while( i != p_transpose.end() ) 00412 { // for each column that this row uses 00413 c_used.begin(i); 00414 ell = c_used.next_element(); 00415 while( ell != c_used.end() ) 00416 { // if this is not the same column, forbid its color 00417 if( ell < j ) 00418 forbidden[ color[ell] ] = true; 00419 ell = c_used.next_element(); 00420 } 00421 i = p_transpose.next_element(); 00422 } 00423 00424 // for each row that this column uses 00425 r_used.begin(j); 00426 i = r_used.next_element(); 00427 while( i != r_used.end() ) 00428 { // For each column that is non-zero for this row 00429 // (the used columns have already been checked above). 00430 not_used.begin(i); 00431 ell = not_used.next_element(); 00432 while( ell != not_used.end() ) 00433 { // if this is not the same column, forbid its color 00434 if( ell < j ) 00435 forbidden[ color[ell] ] = true; 00436 ell = not_used.next_element(); 00437 } 00438 i = r_used.next_element(); 00439 } 00440 00441 // pick the color with smallest index 00442 ell = 0; 00443 while( forbidden[ell] ) 00444 { ell++; 00445 CPPAD_ASSERT_UNKNOWN( ell <= j ); 00446 } 00447 color[j] = ell; 00448 } 00449 } 00450 size_t n_color = 1; 00451 for(ell = 0; ell < n; ell++) 00452 n_color = std::max(n_color, color[ell] + 1); 00453 00454 // direction vector for calls to forward 00455 VectorBase dx(n); 00456 00457 // location for return values from forward 00458 VectorBase dy(m); 00459 00460 // initialize the return value 00461 for(k = 0; k < K; k++) 00462 jac[k] = zero; 00463 00464 // loop over colors 00465 size_t n_sweep = 0; 00466 for(ell = 0; ell < n_color; ell++) 00467 { bool any = false; 00468 k = 0; 00469 for(j = 0; j < n; j++) if( color[j] == ell ) 00470 { // find first k such that col[sort_col[k]] has color ell 00471 if( ! any ) 00472 { while( col[sort_col[k]] < j ) 00473 k++; 00474 any = col[sort_col[k]] == j; 00475 } 00476 } 00477 if( any ) 00478 { n_sweep++; 00479 // combine all columns with this color 00480 for(j = 0; j < n; j++) 00481 { dx[j] = zero; 00482 if( color[j] == ell ) 00483 dx[j] = one; 00484 } 00485 // call forward mode for all these columns at once 00486 dy = Forward(1, dx); 00487 00488 // set the corresponding components of the result 00489 for(j = 0; j < n; j++) if( color[j] == ell ) 00490 { // find first index in c for this column 00491 while( col[sort_col[k]] < j ) 00492 k++; 00493 // extract the row results for this column 00494 while( col[sort_col[k]] == j ) 00495 { jac[ sort_col[k] ] = dy[row[sort_col[k]]]; 00496 k++; 00497 } 00498 } 00499 } 00500 } 00501 return n_sweep; 00502 } 00503 /*! 00504 Private helper function for reverse mode cases. 00505 00506 \tparam Base 00507 See \c SparseJacobianForward(x, p, row, col, jac, work). 00508 00509 \tparam VectorBase 00510 See \c SparseJacobianForward(x, p, row, col, jac, work). 00511 00512 \tparam VectorSet 00513 is either \c sparse_pack, \c sparse_set or \c sparse_list. 00514 00515 \param x 00516 See \c SparseJacobianForward(x, p, row, col, jac, work). 00517 00518 \param p 00519 If <code>work.color.size() != 0</code>, then \c p is not used. 00520 Otherwise, it is a 00521 sparsity pattern for the Jacobian of this ADFun<Base> object. 00522 Note that we do not change the values in \c p_transpose, 00523 but is not \c const because we use its iterator facility. 00524 00525 \param jac 00526 See \c SparseJacobianForward(x, p, row, col, jac, work). 00527 00528 \param work 00529 See \c SparseJacobianReverse(x, p, row, col, jac, work). 00530 In addition, \c work.color contains information that is computed by 00531 \c SparseJacobainRev. 00532 If the sparsity pattern, \c row vector, and \c col vectors 00533 are not the same between calls, \c work.clear() must be called 00534 to reinitialize \c work.color. 00535 00536 \return 00537 See \c SparseJacobianReverse(x, p, row, col, jac, work). 00538 */ 00539 template<class Base> 00540 template <class VectorBase, class VectorSet> 00541 size_t ADFun<Base>::SparseJacobianRev( 00542 const VectorBase& x , 00543 VectorSet& p , 00544 VectorBase& jac , 00545 sparse_jacobian_work& work ) 00546 { 00547 using CppAD::vectorBool; 00548 size_t i, j, k, ell; 00549 00550 CppAD::vector<size_t>& row(work.user_row); 00551 CppAD::vector<size_t>& col(work.user_col); 00552 CppAD::vector<size_t>& sort_row(work.sort_row); 00553 CppAD::vector<size_t>& color(work.color); 00554 00555 size_t m = Range(); 00556 size_t n = Domain(); 00557 00558 // some values 00559 const Base zero(0); 00560 const Base one(1); 00561 00562 // check VectorBase is Simple Vector class with Base type elements 00563 CheckSimpleVector<Base, VectorBase>(); 00564 00565 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00566 CPPAD_ASSERT_UNKNOWN (color.size() == m || color.size() == 0 ); 00567 00568 // number of components of Jacobian that are required 00569 size_t K = size_t(jac.size()); 00570 CPPAD_ASSERT_UNKNOWN( row.size() == K+1 ); 00571 CPPAD_ASSERT_UNKNOWN( col.size() == K+1 ); 00572 CPPAD_ASSERT_UNKNOWN( row[K] == m ); 00573 CPPAD_ASSERT_UNKNOWN( col[K] == n ); 00574 00575 // Point at which we are evaluating the Jacobian 00576 Forward(0, x); 00577 00578 if( color.size() == 0 ) 00579 { CPPAD_ASSERT_UNKNOWN( p.n_set() == m ); 00580 CPPAD_ASSERT_UNKNOWN( p.end() == n ); 00581 00582 // rows and columns that are in the returned jacobian 00583 VectorSet r_used, c_used; 00584 r_used.resize(n, m); 00585 c_used.resize(m, n); 00586 k = 0; 00587 while( k < K ) 00588 { CPPAD_ASSERT_UNKNOWN( 00589 row[sort_row[k]] < m && col[sort_row[k]] < n 00590 ); 00591 CPPAD_ASSERT_UNKNOWN( 00592 k == 0 || row[sort_row[k-1]] <= row[sort_row[k]] 00593 ); 00594 CPPAD_ASSERT_KNOWN( 00595 p.is_element(row[sort_row[k]], col[sort_row[k]]) , 00596 "SparseJacobianReverse: " 00597 "an (row, col) pair is not in sparsity pattern." 00598 ); 00599 r_used.add_element(col[sort_row[k]], row[sort_row[k]]); 00600 c_used.add_element(row[sort_row[k]], col[sort_row[k]]); 00601 k++; 00602 } 00603 00604 // given a column index, which rows are non-zero and not used 00605 VectorSet not_used; 00606 not_used.resize(n, m); 00607 for(i = 0; i < m; i++) 00608 { p.begin(i); 00609 j = p.next_element(); 00610 while( j != p.end() ) 00611 { if( ! r_used.is_element(j , i) ) 00612 not_used.add_element(j, i); 00613 j = p.next_element(); 00614 } 00615 } 00616 00617 // initial coloring 00618 color.resize(m); 00619 for(i = 0; i < m; i++) 00620 color[i] = i; 00621 00622 // See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00623 // Graph Coloring in Optimization Revisited by 00624 // Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00625 vectorBool forbidden(m); 00626 for(i = 1; i < m; i++) 00627 { 00628 // initial all colors as ok for this row 00629 // (value of forbidden for ell > i does not matter) 00630 for(ell = 0; ell <= i; ell++) 00631 forbidden[ell] = false; 00632 00633 // ----------------------------------------------------- 00634 // Forbid colors for which this row would destroy results 00635 // for each column that is non-zero for this row 00636 p.begin(i); 00637 j = p.next_element(); 00638 while( j != p.end() ) 00639 { // for each row that this column uses 00640 r_used.begin(j); 00641 ell = r_used.next_element(); 00642 while( ell != r_used.end() ) 00643 { // if this is not the same row, forbid its color 00644 if( ell < i ) 00645 forbidden[ color[ell] ] = true; 00646 ell = r_used.next_element(); 00647 } 00648 j = p.next_element(); 00649 } 00650 00651 00652 // ----------------------------------------------------- 00653 // Forbid colors that would destroy results for this row. 00654 // for each column that this row uses 00655 c_used.begin(i); 00656 j = c_used.next_element(); 00657 while( j != c_used.end() ) 00658 { // For each row that is non-zero for this column 00659 // (the used rows have already been checked above). 00660 not_used.begin(j); 00661 ell = not_used.next_element(); 00662 while( ell != not_used.end() ) 00663 { // if this is not the same row, forbid its color 00664 if( ell < i ) 00665 forbidden[ color[ell] ] = true; 00666 ell = not_used.next_element(); 00667 } 00668 j = c_used.next_element(); 00669 } 00670 00671 // pick the color with smallest index 00672 ell = 0; 00673 while( forbidden[ell] ) 00674 { ell++; 00675 CPPAD_ASSERT_UNKNOWN( ell <= i ); 00676 } 00677 color[i] = ell; 00678 } 00679 } 00680 size_t n_color = 1; 00681 for(ell = 0; ell < m; ell++) 00682 n_color = std::max(n_color, color[ell] + 1); 00683 00684 // weighting vector for calls to reverse 00685 VectorBase w(m); 00686 00687 // location for return values from Reverse 00688 VectorBase dw(n); 00689 00690 // initialize the return value 00691 for(k = 0; k < K; k++) 00692 jac[k] = zero; 00693 00694 // loop over colors 00695 size_t n_sweep = 0; 00696 for(ell = 0; ell < n_color; ell++) 00697 { bool any = false; 00698 k = 0; 00699 for(i = 0; i < m; i++) if( color[i] == ell ) 00700 { // find first k such that row[sort_row[k]] has color ell 00701 if( ! any ) 00702 { while( row[sort_row[k]] < i ) 00703 k++; 00704 any = row[sort_row[k]] == i; 00705 } 00706 } 00707 if( any ) 00708 { n_sweep++; 00709 // combine all the rows with this color 00710 for(i = 0; i < m; i++) 00711 { w[i] = zero; 00712 if( color[i] == ell ) 00713 w[i] = one; 00714 } 00715 // call reverse mode for all these rows at once 00716 dw = Reverse(1, w); 00717 00718 // set the corresponding components of the result 00719 for(i = 0; i < m; i++) if( color[i] == ell ) 00720 { // find first index in r for this row 00721 while( row[sort_row[k]] < i ) 00722 k++; 00723 // extract the row results for this row 00724 while( row[sort_row[k]] == i ) 00725 { jac[ sort_row[k] ] = dw[col[sort_row[k]]]; 00726 k++; 00727 } 00728 } 00729 } 00730 } 00731 return n_sweep; 00732 } 00733 // ========================================================================== 00734 // Public Member functions 00735 // ========================================================================== 00736 /*! 00737 Compute user specified subset of a sparse Jacobian using forward mode. 00738 00739 The C++ source code corresponding to this operation is 00740 \verbatim 00741 SparceJacobianForward(x, p, row, col, jac, work) 00742 \endverbatim 00743 00744 \tparam Base 00745 is the base type for the recording that is stored in this 00746 ADFun<Base object. 00747 00748 \tparam VectorBase 00749 is a simple vector class with elements of type \a Base. 00750 00751 \tparam VectorSet 00752 is a simple vector class with elements of type 00753 \c bool or \c std::set<size_t>. 00754 00755 \tparam VectorSize 00756 is a simple vector class with elements of type \c size_t. 00757 00758 \param x 00759 is a vector specifing the point at which to compute the Jacobian. 00760 00761 \param p 00762 is the sparsity pattern for the Jacobian that we are calculating. 00763 00764 \param row 00765 is the vector of row indices for the returned Jacobian values. 00766 00767 \param col 00768 is the vector of columns indices for the returned Jacobian values. 00769 It must have the same size are r. 00770 00771 \param jac 00772 is the vector of Jacobian values. 00773 It must have the same size are r. 00774 The return value <code>jac[k]</code> is the partial of the 00775 <code>row[k]</code> component of the function with respect 00776 the the <code>col[k]</code> of its argument. 00777 00778 \param work 00779 \c work contains information that depends on the function object, 00780 sparsity pattern, \c row vector, and \c col vector. 00781 These \c work values are computed by \c SparseJacobianForward. 00782 If they are not the same between calls, \c work.clear() must be called 00783 to reinitialize them. 00784 To be more specific: 00785 \c user_row is a copy of \c row with the value \c m at the end 00786 (\c m is the range dimension), 00787 \c user_col is a copy of \c col with the value \c n at the end 00788 (\c n is the domain dimension), 00789 \c sort_col is the index vector that sorts the \c user_col 00790 vector with the extra value \c K at the end 00791 (\c K is the size of \c jac). 00792 00793 \return 00794 Is the number of first order forward sweeps used to compute the 00795 requested Jacobian values. The total work, not counting the zero order 00796 forward sweep, or the time to combine computations, is proportional to this 00797 return value. 00798 */ 00799 template<class Base> 00800 template <class VectorBase, class VectorSet, class VectorSize> 00801 size_t ADFun<Base>::SparseJacobianForward( 00802 const VectorBase& x , 00803 const VectorSet& p , 00804 const VectorSize& row , 00805 const VectorSize& col , 00806 VectorBase& jac , 00807 sparse_jacobian_work& work ) 00808 { 00809 size_t n = Domain(); 00810 size_t m = Range(); 00811 size_t k, K = jac.size(); 00812 if( work.user_row.size() == 0 ) 00813 { // create version of (row, col, k) sorted by column value 00814 work.user_col.resize(K+1); 00815 work.user_row.resize(K+1); 00816 work.sort_col.resize(K+1); 00817 00818 // put sorted indices in user_row and user_col 00819 for(k = 0; k < K; k++) 00820 { work.user_row[k] = row[k]; 00821 work.user_col[k] = col[k]; 00822 } 00823 work.user_row[K] = m; 00824 work.user_col[K] = n; 00825 00826 // put sorting indices in sort_col 00827 index_sort(work.user_col, work.sort_col); 00828 } 00829 # ifndef NDEBUG 00830 CPPAD_ASSERT_KNOWN( 00831 size_t(x.size()) == n , 00832 "SparseJacobianForward: size of x not equal domain dimension for f." 00833 ); 00834 CPPAD_ASSERT_KNOWN( 00835 size_t(row.size()) == K && size_t(col.size()) == K , 00836 "SparseJacobianForward: either r or c does not have " 00837 "the same size as jac." 00838 ); 00839 CPPAD_ASSERT_KNOWN( 00840 work.user_row.size() == K+1 && 00841 work.user_col.size() == K+1 && 00842 work.sort_col.size() == K+1 , 00843 "SparseJacobianForward: invalid value in work." 00844 ); 00845 CPPAD_ASSERT_KNOWN( 00846 work.color.size() == 0 || work.color.size() == n, 00847 "SparseJacobianForward: invalid value in work." 00848 ); 00849 for(k = 0; k < K; k++) 00850 { CPPAD_ASSERT_KNOWN( 00851 row[k] < m, 00852 "SparseJacobianForward: invalid value in r." 00853 ); 00854 CPPAD_ASSERT_KNOWN( 00855 col[k] < n, 00856 "SparseJacobianForward: invalid value in c." 00857 ); 00858 CPPAD_ASSERT_KNOWN( 00859 work.sort_col[k] < K, 00860 "SparseJacobianForward: invalid value in work." 00861 ); 00862 CPPAD_ASSERT_KNOWN( 00863 work.user_row[k] == row[k], 00864 "SparseJacobianForward: invalid value in work." 00865 ); 00866 CPPAD_ASSERT_KNOWN( 00867 work.user_col[k] == col[k], 00868 "SparseJacobianForward: invalid value in work." 00869 ); 00870 } 00871 if( work.color.size() != 0 ) 00872 for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN( 00873 work.color[j] < n, 00874 "SparseJacobianForward: invalid value in work." 00875 ); 00876 # endif 00877 00878 typedef typename VectorSet::value_type Set_type; 00879 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00880 Pattern_type s_transpose; 00881 if( work.color.size() == 0 ) 00882 { bool transpose = true; 00883 sparsity_user2internal(s_transpose, p, m, n, transpose); 00884 } 00885 size_t n_sweep = SparseJacobianFor(x, s_transpose, jac, work); 00886 return n_sweep; 00887 } 00888 /*! 00889 Compute user specified subset of a sparse Jacobian using forward mode. 00890 00891 The C++ source code corresponding to this operation is 00892 \verbatim 00893 SparceJacobianReverse(x, p, row, col, jac, work) 00894 \endverbatim 00895 00896 \tparam Base 00897 See \c SparseJacobianForward(x, p, row, col, jac, work). 00898 00899 \tparam VectorBase 00900 See \c SparseJacobianForward(x, p, row, col, jac, work). 00901 00902 \tparam VectorSet 00903 See \c SparseJacobianForward(x, p, row, col, jac, work). 00904 00905 \tparam VectorSize 00906 See \c SparseJacobianForward(x, p, row, col, jac, work). 00907 00908 \param x 00909 See \c SparseJacobianForward(x, p, row, col, jac, work). 00910 00911 \param p 00912 See \c SparseJacobianForward(x, p, row, col, jac, work). 00913 00914 \param row 00915 See \c SparseJacobianForward(x, p, row, col, jac, work). 00916 00917 \param col 00918 See \c SparseJacobianForward(x, p, row, col, jac, work). 00919 00920 \param jac 00921 See \c SparseJacobianForward(x, p, row, col, jac, work). 00922 00923 \param work 00924 \c work contains information that depends on the function object, 00925 sparsity pattern, \c row vector, and \c col vector. 00926 These \c work values are computed by \c SparseJacobianReverse. 00927 If they are not the same between calls, \c work.clear() must be called 00928 to reinitialize them. 00929 To be more specific: 00930 \c user_row is a copy of \c row with the value \c m at the end 00931 (\c m is the range dimension), 00932 \c user_col is a copy of \c col with the value \c n at the end 00933 (\c n is the domain dimension), 00934 \c sort_row is the index vector that sorts the \c user_row 00935 vector with the extra value \c K at the end 00936 (\c K is the size of \c jac). 00937 00938 \return 00939 Is the number of first order reverse sweeps used to compute the 00940 reverse Jacobian values. The total work, not counting the zero order 00941 forward sweep, or the time to combine computations, is proportional to this 00942 return value. 00943 */ 00944 template<class Base> 00945 template <class VectorBase, class VectorSet, class VectorSize> 00946 size_t ADFun<Base>::SparseJacobianReverse( 00947 const VectorBase& x , 00948 const VectorSet& p , 00949 const VectorSize& row , 00950 const VectorSize& col , 00951 VectorBase& jac , 00952 sparse_jacobian_work& work ) 00953 { 00954 size_t m = Range(); 00955 size_t n = Domain(); 00956 size_t k, K = jac.size(); 00957 if( work.user_row.size() == 0 ) 00958 { // create version of (row, col, k) sorted by row row value 00959 work.user_col.resize(K+1); 00960 work.user_row.resize(K+1); 00961 work.sort_row.resize(K+1); 00962 00963 // put sorted indices in user_row and user_col 00964 for(k = 0; k < K; k++) 00965 { work.user_row[k] = row[k]; 00966 work.user_col[k] = col[k]; 00967 } 00968 work.user_row[K] = m; 00969 work.user_col[K] = n; 00970 00971 // put sorting indices in sort_row 00972 index_sort(work.user_row, work.sort_row); 00973 } 00974 # ifndef NDEBUG 00975 CPPAD_ASSERT_KNOWN( 00976 size_t(x.size()) == n , 00977 "SparseJacobianReverse: size of x not equal domain dimension for f." 00978 ); 00979 CPPAD_ASSERT_KNOWN( 00980 size_t(row.size()) == K && size_t(col.size()) == K , 00981 "SparseJacobianReverse: either r or c does not have " 00982 "the same size as jac." 00983 ); 00984 CPPAD_ASSERT_KNOWN( 00985 work.user_row.size() == K+1 && 00986 work.user_col.size() == K+1 && 00987 work.sort_row.size() == K+1 , 00988 "SparseJacobianReverse: invalid value in work." 00989 ); 00990 CPPAD_ASSERT_KNOWN( 00991 work.color.size() == 0 || work.color.size() == m, 00992 "SparseJacobianReverse: invalid value in work." 00993 ); 00994 for(k = 0; k < K; k++) 00995 { CPPAD_ASSERT_KNOWN( 00996 row[k] < m, 00997 "SparseJacobianReverse: invalid value in r." 00998 ); 00999 CPPAD_ASSERT_KNOWN( 01000 col[k] < n, 01001 "SparseJacobianReverse: invalid value in c." 01002 ); 01003 CPPAD_ASSERT_KNOWN( 01004 work.sort_row[k] < K, 01005 "SparseJacobianReverse: invalid value in work." 01006 ); 01007 CPPAD_ASSERT_KNOWN( 01008 work.user_row[k] == row[k], 01009 "SparseJacobianReverse: invalid value in work." 01010 ); 01011 CPPAD_ASSERT_KNOWN( 01012 work.user_col[k] == col[k], 01013 "SparseJacobianReverse: invalid value in work." 01014 ); 01015 } 01016 if( work.color.size() != 0 ) 01017 for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN( 01018 work.color[i] < m, 01019 "SparseJacobianReverse: invalid value in work." 01020 ); 01021 # endif 01022 01023 typedef typename VectorSet::value_type Set_type; 01024 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 01025 Pattern_type s; 01026 if( work.color.size() == 0 ) 01027 { bool transpose = false; 01028 sparsity_user2internal(s, p, m, n, transpose); 01029 } 01030 size_t n_sweep = SparseJacobianRev(x, s, jac, work); 01031 return n_sweep; 01032 } 01033 /*! 01034 Compute a sparse Jacobian. 01035 01036 The C++ source code corresponding to this operation is 01037 \verbatim 01038 jac = SparseJacobian(x, p) 01039 \endverbatim 01040 01041 \tparam Base 01042 is the base type for the recording that is stored in this 01043 ADFun<Base object. 01044 01045 \tparam VectorBase 01046 is a simple vector class with elements of type \a Base. 01047 01048 \tparam VectorSet 01049 is a simple vector class with elements of type 01050 \c bool or \c std::set<size_t>. 01051 01052 \param x 01053 is a vector specifing the point at which to compute the Jacobian. 01054 01055 \param p 01056 is the sparsity pattern for the Jacobian that we are calculating. 01057 01058 \return 01059 Will be a vector if size \c m * n containing the Jacobian at the 01060 specified point (in row major order). 01061 */ 01062 template <class Base> 01063 template <class VectorBase, class VectorSet> 01064 VectorBase ADFun<Base>::SparseJacobian( 01065 const VectorBase& x, const VectorSet& p 01066 ) 01067 { size_t i, j, k; 01068 01069 size_t m = Range(); 01070 size_t n = Domain(); 01071 VectorBase jac(m * n); 01072 01073 CPPAD_ASSERT_KNOWN( 01074 size_t(x.size()) == n, 01075 "SparseJacobian: size of x not equal domain size for f." 01076 ); 01077 CheckSimpleVector<Base, VectorBase>(); 01078 01079 typedef typename VectorSet::value_type Set_type; 01080 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 01081 01082 // initialize the return value as zero 01083 Base zero(0); 01084 for(i = 0; i < m; i++) 01085 for(j = 0; j < n; j++) 01086 jac[i * n + j] = zero; 01087 01088 sparse_jacobian_work work; 01089 if( n <= m ) 01090 { CppAD::vector<size_t>& row(work.user_row); 01091 CppAD::vector<size_t>& col(work.user_col); 01092 CppAD::vector<size_t>& sort_col(work.sort_col); 01093 01094 // forward mode, columns are sorted 01095 assert( row.size() == 0 ); 01096 assert( col.size() == 0 ); 01097 assert( sort_col.size() == 0 ); 01098 01099 // need an internal copy of sparsity pattern 01100 Pattern_type s_transpose; 01101 bool transpose = true; 01102 sparsity_user2internal(s_transpose, p, m, n, transpose); 01103 01104 k = 0; 01105 for(j = 0; j < n; j++) 01106 { s_transpose.begin(j); 01107 i = s_transpose.next_element(); 01108 while( i != s_transpose.end() ) 01109 { row.push_back(i); 01110 col.push_back(j); 01111 sort_col.push_back(k); 01112 k++; 01113 i = s_transpose.next_element(); 01114 } 01115 } 01116 size_t K = k; 01117 VectorBase J(K); 01118 row.push_back(m); 01119 col.push_back(n); 01120 sort_col.push_back(K); 01121 01122 // now we have folded this into the following case 01123 SparseJacobianFor(x, s_transpose, J, work); 01124 01125 // now set the non-zero return values 01126 for(k = 0; k < K; k++) 01127 jac[row[sort_col[k]] * n + col[sort_col[k]]] = J[k]; 01128 } 01129 else 01130 { CppAD::vector<size_t>& row(work.user_row); 01131 CppAD::vector<size_t>& col(work.user_col); 01132 CppAD::vector<size_t>& sort_row(work.sort_row); 01133 01134 // reverse mode, rows are sorted 01135 assert( row.size() == 0 ); 01136 assert( col.size() == 0 ); 01137 assert( sort_row.size() == 0 ); 01138 01139 // need an internal copy of sparsity pattern 01140 Pattern_type s; 01141 bool transpose = false; 01142 sparsity_user2internal(s, p, m, n, transpose); 01143 01144 k = 0; 01145 for(i = 0; i < m; i++) 01146 { s.begin(i); 01147 j = s.next_element(); 01148 while( j != s.end() ) 01149 { row.push_back(i); 01150 col.push_back(j); 01151 sort_row.push_back(k); 01152 k++; 01153 j = s.next_element(); 01154 } 01155 } 01156 size_t K = k; 01157 VectorBase J(K); 01158 row.push_back(m); 01159 col.push_back(n); 01160 sort_row.push_back(K); 01161 01162 // now we have folded this into the following case 01163 SparseJacobianRev(x, s, J, work); 01164 01165 // now set the non-zero return values 01166 for(k = 0; k < K; k++) 01167 jac[row[sort_row[k]] * n + col[sort_row[k]]] = J[k]; 01168 } 01169 01170 return jac; 01171 } 01172 01173 /*! 01174 Compute a sparse Jacobian. 01175 01176 The C++ source code corresponding to this operation is 01177 \verbatim 01178 jac = SparseJacobian(x) 01179 \endverbatim 01180 01181 \tparam Base 01182 is the base type for the recording that is stored in this 01183 ADFun<Base object. 01184 01185 \tparam VectorBase 01186 is a simple vector class with elements of the \a Base. 01187 01188 \param x 01189 is a vector specifing the point at which to compute the Jacobian. 01190 01191 \return 01192 Will be a vector of size \c m * n containing the Jacobian at the 01193 specified point (in row major order). 01194 */ 01195 template <class Base> 01196 template <class VectorBase> 01197 VectorBase ADFun<Base>::SparseJacobian( const VectorBase& x ) 01198 { typedef CppAD::vectorBool VectorBool; 01199 01200 size_t m = Range(); 01201 size_t n = Domain(); 01202 01203 // sparsity pattern for Jacobian 01204 VectorBool p(m * n); 01205 01206 if( n <= m ) 01207 { size_t j, k; 01208 01209 // use forward mode 01210 VectorBool r(n * n); 01211 for(j = 0; j < n; j++) 01212 { for(k = 0; k < n; k++) 01213 r[j * n + k] = false; 01214 r[j * n + j] = true; 01215 } 01216 p = ForSparseJac(n, r); 01217 } 01218 else 01219 { size_t i, k; 01220 01221 // use reverse mode 01222 VectorBool s(m * m); 01223 for(i = 0; i < m; i++) 01224 { for(k = 0; k < m; k++) 01225 s[i * m + k] = false; 01226 s[i * m + i] = true; 01227 } 01228 p = RevSparseJac(m, s); 01229 } 01230 return SparseJacobian(x, p); 01231 } 01232 01233 01234 /*! \} */ 01235 CPPAD_END_NAMESPACE 01236 # endif