CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_SPARSE_HESSIAN_INCLUDED 00003 # define CPPAD_SPARSE_HESSIAN_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 /* 00017 $begin sparse_hessian$$ 00018 $spell 00019 recomputed 00020 CppAD 00021 valarray 00022 std 00023 Bool 00024 hes 00025 const 00026 Taylor 00027 $$ 00028 00029 $section Sparse Hessian: Easy Driver$$ 00030 $index SparseHessian$$ 00031 $index hessian, sparse$$ 00032 00033 $head Syntax$$ 00034 $icode%hes% = %f%.SparseHessian(%x%, %w%) 00035 %hes% = %f%.SparseHessian(%x%, %w%, %p%) 00036 %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%) 00037 %$$ 00038 00039 $head Purpose$$ 00040 We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size, 00041 and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$. 00042 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the 00043 $cref/AD function/glossary/AD Function/$$ 00044 corresponding to $icode f$$. 00045 The syntax above sets $icode hes$$ to the Hessian 00046 $latex \[ 00047 H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x) 00048 \] $$ 00049 This routine takes advantage of the sparsity of the Hessian 00050 in order to reduce the amount of computation necessary. 00051 If $icode row$$ and $icode col$$ are present, it also takes 00052 advantage of the reduced set of elements of the Hessian that 00053 need to be computed. 00054 One can use speed tests (e.g. $cref speed_test$$) 00055 to verify that results are computed faster 00056 than when using the routine $cref Hessian$$. 00057 00058 $head f$$ 00059 The object $icode f$$ has prototype 00060 $codei% 00061 ADFun<%Base%> %f% 00062 %$$ 00063 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$ 00064 (see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below). 00065 00066 $head x$$ 00067 The argument $icode x$$ has prototype 00068 $codei% 00069 const %VectorBase%& %x% 00070 %$$ 00071 (see $cref/VectorBase/sparse_hessian/VectorBase/$$ below) 00072 and its size 00073 must be equal to $icode n$$, the dimension of the 00074 $cref/domain/seq_property/Domain/$$ space for $icode f$$. 00075 It specifies 00076 that point at which to evaluate the Hessian. 00077 00078 $head w$$ 00079 The argument $icode w$$ has prototype 00080 $codei% 00081 const %VectorBase%& %w% 00082 %$$ 00083 and size $latex m$$. 00084 It specifies the value of $latex w_i$$ in the expression 00085 for $icode hes$$. 00086 The more components of $latex w$$ that are identically zero, 00087 the more sparse the resulting Hessian may be (and hence the more efficient 00088 the calculation of $icode hes$$ may be). 00089 00090 $head p$$ 00091 The argument $icode p$$ is optional and has prototype 00092 $codei% 00093 const %VectorSet%& %p% 00094 %$$ 00095 (see $cref/VectorSet/sparse_hessian/VectorSet/$$ below) 00096 If it has elements of type $code bool$$, 00097 its size is $latex n * n$$. 00098 If it has elements of type $code std::set<size_t>$$, 00099 its size is $latex n$$ and all its set elements are between 00100 zero and $latex n - 1$$. 00101 It specifies a 00102 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 00103 for the Hessian $latex H(x)$$. 00104 $pre 00105 00106 $$ 00107 If this sparsity pattern does not change between calls to 00108 $codei SparseHessian$$, it should be faster to calculate $icode p$$ once and 00109 pass this argument to $codei SparseHessian$$. 00110 In addition, 00111 if you specify $icode p$$, CppAD will use the same 00112 type of sparsity representation 00113 (vectors of $code bool$$ or vectors of $code std::set<size_t>$$) 00114 for its internal calculations. 00115 Otherwise, the representation 00116 for the internal calculations is unspecified. 00117 00118 $head row, col$$ 00119 The arguments $icode row$$ and $icode col$$ are optional and have prototype 00120 $codei% 00121 const %VectorSize%& %row% 00122 const %VectorSize%& %col% 00123 %$$ 00124 (see $cref/VectorSize/sparse_hessian/VectorSize/$$ below). 00125 They specify which rows and columns of $latex H (x)$$ are 00126 returned and in what order. 00127 We use $latex K$$ to denote the value $icode%hes%.size()%$$ 00128 which must also equal the size of $icode row$$ and $icode col$$. 00129 Furthermore, 00130 for $latex k = 0 , \ldots , K-1$$, it must hold that 00131 $latex row[k] < n$$ and $latex col[k] < n$$. 00132 In addition, 00133 all of the $latex (row[k], col[k])$$ pairs must correspond to a true value 00134 in the sparsity pattern $icode p$$. 00135 00136 $head hes$$ 00137 The result $icode hes$$ has prototype 00138 $codei% 00139 %VectorBase% %hes% 00140 %$$ 00141 In the case where $icode row$$ and $icode col$$ are not present, 00142 the size of $icode hes$$ is $latex n * n$$ and 00143 its size is $latex n * n$$. 00144 In this case, for $latex i = 0 , \ldots , n - 1 $$ 00145 and $latex ell = 0 , \ldots , n - 1$$ 00146 $latex \[ 00147 hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x ) 00148 \] $$ 00149 $pre 00150 00151 $$ 00152 In the case where the arguments $icode row$$ and $icode col$$ are present, 00153 we use $latex K$$ to denote the size of $icode hes$$. 00154 The input value of its elements does not matter. 00155 Upon return, for $latex k = 0 , \ldots , K - 1$$, 00156 $latex \[ 00157 hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x) 00158 \; , \; 00159 \; {\rm where} \; 00160 j = row[k] 00161 \; {\rm and } \; 00162 \ell = col[k] 00163 \] $$ 00164 00165 $head work$$ 00166 If this argument is present, it has prototype 00167 $codei% 00168 sparse_hessian_work& %work% 00169 %$$ 00170 This object can only be used with the routines $code SparseHessian$$. 00171 During its the first use, information is stored in $icode work$$. 00172 This is used to reduce the work done by future calls to $code SparseHessian$$ 00173 with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$. 00174 If a future call is make where any of these values have changed, 00175 you must first call $icode%work%.clear()%$$ 00176 to inform CppAD that this information needs to be recomputed. 00177 00178 $head n_sweep$$ 00179 The return value $icode n_sweep$$ has prototype 00180 $codei% 00181 size_t %n_sweep% 00182 %$$ 00183 It is the number of first order forward sweeps 00184 used to compute the requested Hessian values. 00185 Each first forward sweep is followed by a second order reverse sweep 00186 so it is also the number of reverse sweeps. 00187 This is proportional to the total work that $code SparseHessian$$ does, 00188 not counting the zero order forward sweep, 00189 or the work to combine multiple columns into a single 00190 forward-reverse sweep pair. 00191 00192 $head VectorBase$$ 00193 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with 00194 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00195 $icode Base$$. 00196 The routine $cref CheckSimpleVector$$ will generate an error message 00197 if this is not the case. 00198 00199 $head VectorSet$$ 00200 The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with 00201 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00202 $code bool$$ or $code std::set<size_t>$$; 00203 see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion 00204 of the difference. 00205 The routine $cref CheckSimpleVector$$ will generate an error message 00206 if this is not the case. 00207 00208 $subhead Restrictions$$ 00209 If $icode VectorSet$$ has elements of $code std::set<size_t>$$, 00210 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 00211 corresponding set. 00212 According to section 26.3.2.3 of the 1998 C++ standard, 00213 $code std::valarray< std::set<size_t> >$$ does not satisfy 00214 this condition. 00215 00216 $head VectorSize$$ 00217 The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with 00218 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00219 $code size_t$$. 00220 The routine $cref CheckSimpleVector$$ will generate an error message 00221 if this is not the case. 00222 00223 $head Uses Forward$$ 00224 After each call to $cref Forward$$, 00225 the object $icode f$$ contains the corresponding 00226 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 00227 After a call to any of the sparse Hessian routines, 00228 the zero order Taylor coefficients correspond to 00229 $icode%f%.Forward(0, %x%)%$$ 00230 and the other coefficients are unspecified. 00231 00232 $head Example$$ 00233 $children% 00234 example/sparse_hessian.cpp 00235 %$$ 00236 The routine 00237 $cref sparse_hessian.cpp$$ 00238 is examples and tests of $code sparse_hessian$$. 00239 It return $code true$$, if it succeeds and $code false$$ otherwise. 00240 00241 00242 $end 00243 ----------------------------------------------------------------------------- 00244 */ 00245 # include <cppad/local/std_set.hpp> 00246 00247 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00248 /*! 00249 \file sparse_hessian.hpp 00250 Sparse Hessian driver routine and helper functions. 00251 */ 00252 // =========================================================================== 00253 /*! 00254 class used by SparseHessian to hold information 00255 so it does not need to be recomputed. 00256 */ 00257 class sparse_hessian_work { 00258 public: 00259 /// indices that sort the user row and col arrays by color 00260 CppAD::vector<size_t> order; 00261 /// results of the coloring algorithm 00262 CppAD::vector<size_t> color; 00263 /// inform CppAD that this information needs to be recomputed 00264 void clear(void) 00265 { order.clear(); 00266 color.clear(); 00267 } 00268 }; 00269 // =========================================================================== 00270 /*! 00271 Private helper function that does computation for all Sparse Hessian cases. 00272 00273 \tparam Base 00274 is the base type for the recording that is stored in this ADFun<Base object. 00275 00276 \tparam VectorBase 00277 is a simple vector class with elements of type \a Base. 00278 00279 \tparam VectorSet 00280 is a simple vector class with elements of type 00281 \c bool or \c std::set<size_t>. 00282 00283 \tparam VectorSize 00284 is \c sparse_pack, \c sparse_set or \c sparse_list. 00285 00286 \param x [in] 00287 is a vector specifing the point at which to compute the Hessian. 00288 00289 \param w [in] 00290 is the weighting vector that defines a scalar valued function by 00291 a weighted sum of the components of the vector valued function 00292 $latex F(x)$$. 00293 00294 \param sparsity [in] 00295 is the sparsity pattern for the Hessian that we are calculating. 00296 00297 \param row [in] 00298 is the vector of row indices for the returned Hessian values. 00299 00300 \param col [in] 00301 is the vector of columns indices for the returned Hessian values. 00302 It must have the same size are r. 00303 00304 \param hes [out] 00305 is the vector of Hessian values. 00306 It must have the same size are r. 00307 The return value <code>hes[k]</code> is the second partial of 00308 \f$ w^{\rm T} F(x)\f$ with respect to the 00309 <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. 00310 00311 \param work 00312 This structure contains information that is computed by \c SparseHessianCompute. 00313 If the sparsity pattern, \c row vector, or \c col vectors 00314 are not the same between calls to \c SparseHessianCompute, 00315 \c work.clear() must be called to reinitialize \c work. 00316 00317 \return 00318 Is the number of first order forward sweeps used to compute the 00319 requested Hessian values. 00320 (This is also equal to the number of second order reverse sweeps.) 00321 The total work, not counting the zero order 00322 forward sweep, or the time to combine computations, is proportional to this 00323 return value. 00324 */ 00325 template<class Base> 00326 template <class VectorBase, class VectorSet, class VectorSize> 00327 size_t ADFun<Base>::SparseHessianCompute( 00328 const VectorBase& x , 00329 const VectorBase& w , 00330 VectorSet& sparsity , 00331 const VectorSize& row , 00332 const VectorSize& col , 00333 VectorBase& hes , 00334 sparse_hessian_work& work ) 00335 { 00336 using CppAD::vectorBool; 00337 size_t i, k, ell; 00338 00339 CppAD::vector<size_t>& color(work.color); 00340 CppAD::vector<size_t>& order(work.order); 00341 00342 size_t n = Domain(); 00343 00344 // some values 00345 const Base zero(0); 00346 const Base one(1); 00347 00348 // check VectorBase is Simple Vector class with Base type elements 00349 CheckSimpleVector<Base, VectorBase>(); 00350 00351 CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n ); 00352 CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n ); 00353 00354 // number of components of Hessian that are required 00355 size_t K = hes.size(); 00356 CPPAD_ASSERT_UNKNOWN( row.size() == K ); 00357 CPPAD_ASSERT_UNKNOWN( col.size() == K ); 00358 00359 // Point at which we are evaluating the Hessian 00360 Forward(0, x); 00361 00362 // check for case where nothing (except Forward above) to do 00363 if( K == 0 ) 00364 return 0; 00365 00366 // Rows of the Hessian (i below) correspond to the forward mode index 00367 // and columns (j below) correspond to the reverse mode index. 00368 if( color.size() == 0 ) 00369 { 00370 CPPAD_ASSERT_UNKNOWN( sparsity.n_set() == n ); 00371 CPPAD_ASSERT_UNKNOWN( sparsity.end() == n ); 00372 00373 // execute coloring algorithm 00374 color.resize(n); 00375 color_general_cppad(sparsity, row, col, color); 00376 00377 // put sorting indices in color order 00378 VectorSize key(K); 00379 order.resize(K); 00380 for(k = 0; k < K; k++) 00381 key[k] = color[ row[k] ]; 00382 index_sort(key, order); 00383 00384 } 00385 size_t n_color = 1; 00386 for(ell = 0; ell < n; ell++) if( color[ell] < n ) 00387 n_color = std::max(n_color, color[ell] + 1); 00388 00389 // direction vector for calls to forward (rows of the Hessian) 00390 VectorBase u(n); 00391 00392 // location for return values from reverse (columns of the Hessian) 00393 VectorBase ddw(2 * n); 00394 00395 // initialize the return value 00396 for(k = 0; k < K; k++) 00397 hes[k] = zero; 00398 00399 // loop over colors 00400 k = 0; 00401 for(ell = 0; ell < n_color; ell++) 00402 { CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell ); 00403 00404 // combine all rows with this color 00405 for(i = 0; i < n; i++) 00406 { u[i] = zero; 00407 if( color[i] == ell ) 00408 u[i] = one; 00409 } 00410 // call forward mode for all these rows at once 00411 Forward(1, u); 00412 00413 // evaluate derivative of w^T * F'(x) * u 00414 ddw = Reverse(2, w); 00415 00416 // set the corresponding components of the result 00417 while( k < K && color[ row[ order[k] ] ] == ell ) 00418 { hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ]; 00419 k++; 00420 } 00421 } 00422 return n_color; 00423 } 00424 // =========================================================================== 00425 // Public Member Functions 00426 // =========================================================================== 00427 /*! 00428 Compute user specified subset of a sparse Hessian. 00429 00430 The C++ source code corresponding to this operation is 00431 \verbatim 00432 SparceHessian(x, w, p, row, col, hes, work) 00433 \endverbatim 00434 00435 \tparam Base 00436 is the base type for the recording that is stored in this ADFun<Base object. 00437 00438 \tparam VectorBase 00439 is a simple vector class with elements of type \a Base. 00440 00441 \tparam VectorSet 00442 is a simple vector class with elements of type 00443 \c bool or \c std::set<size_t>. 00444 00445 \tparam VectorSize 00446 is a simple vector class with elements of type \c size_t. 00447 00448 \param x [in] 00449 is a vector specifing the point at which to compute the Hessian. 00450 00451 \param w [in] 00452 is the weighting vector that defines a scalar valued function by 00453 a weighted sum of the components of the vector valued function 00454 $latex F(x)$$. 00455 00456 \param p [in] 00457 is the sparsity pattern for the Hessian that we are calculating. 00458 00459 \param row [in] 00460 is the vector of row indices for the returned Hessian values. 00461 00462 \param col [in] 00463 is the vector of columns indices for the returned Hessian values. 00464 It must have the same size are r. 00465 00466 \param hes [out] 00467 is the vector of Hessian values. 00468 It must have the same size are r. 00469 The return value <code>hes[k]</code> is the second partial of 00470 \f$ w^{\rm T} F(x)\f$ with respect to the 00471 <code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$. 00472 00473 \param work 00474 This structure contains information that is computed by \c SparseHessianCompute. 00475 If the sparsity pattern, \c row vector, or \c col vectors 00476 are not the same between calls to \c SparseHessian, 00477 \c work.clear() must be called to reinitialize \c work. 00478 00479 \return 00480 Is the number of first order forward sweeps used to compute the 00481 requested Hessian values. 00482 (This is also equal to the number of second order reverse sweeps.) 00483 The total work, not counting the zero order 00484 forward sweep, or the time to combine computations, is proportional to this 00485 return value. 00486 */ 00487 template<class Base> 00488 template <class VectorBase, class VectorSet, class VectorSize> 00489 size_t ADFun<Base>::SparseHessian( 00490 const VectorBase& x , 00491 const VectorBase& w , 00492 const VectorSet& p , 00493 const VectorSize& row , 00494 const VectorSize& col , 00495 VectorBase& hes , 00496 sparse_hessian_work& work ) 00497 { 00498 size_t n = Domain(); 00499 # ifndef NDEBUG 00500 size_t k, K = hes.size(); 00501 CPPAD_ASSERT_KNOWN( 00502 size_t(x.size()) == n , 00503 "SparseHessian: size of x not equal domain dimension for f." 00504 ); 00505 CPPAD_ASSERT_KNOWN( 00506 size_t(row.size()) == K && size_t(col.size()) == K , 00507 "SparseHessian: either r or c does not have the same size as ehs." 00508 ); 00509 CPPAD_ASSERT_KNOWN( 00510 work.color.size() == 0 || work.color.size() == n, 00511 "SparseHessian: invalid value in work." 00512 ); 00513 for(k = 0; k < K; k++) 00514 { CPPAD_ASSERT_KNOWN( 00515 row[k] < n, 00516 "SparseHessian: invalid value in r." 00517 ); 00518 CPPAD_ASSERT_KNOWN( 00519 col[k] < n, 00520 "SparseHessian: invalid value in c." 00521 ); 00522 } 00523 if( work.color.size() != 0 ) 00524 for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN( 00525 work.color[j] <= n, 00526 "SparseHessian: invalid value in work." 00527 ); 00528 # endif 00529 typedef typename VectorSet::value_type Set_type; 00530 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00531 Pattern_type s; 00532 if( work.color.size() == 0 ) 00533 { bool transpose = false; 00534 sparsity_user2internal(s, p, n, n, transpose); 00535 } 00536 size_t n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work); 00537 return n_sweep; 00538 } 00539 /*! 00540 Compute a sparse Hessian. 00541 00542 The C++ source code coresponding to this operation is 00543 \verbatim 00544 hes = SparseHessian(x, w, p) 00545 \endverbatim 00546 00547 00548 \tparam Base 00549 is the base type for the recording that is stored in this 00550 ADFun<Base object. 00551 00552 \tparam VectorBase 00553 is a simple vector class with elements of the \a Base. 00554 00555 \tparam VectorSet 00556 is a simple vector class with elements of type 00557 \c bool or \c std::set<size_t>. 00558 00559 \param x [in] 00560 is a vector specifing the point at which to compute the Hessian. 00561 00562 \param w [in] 00563 The Hessian is computed for a weighted sum of the components 00564 of the function corresponding to this ADFun<Base> object. 00565 The argument \a w specifies the weights for each component. 00566 It must have size equal to the range dimension for this ADFun<Base> object. 00567 00568 \param p [in] 00569 is a sparsity pattern for the Hessian. 00570 00571 \return 00572 Will be a vector of size \c n * n containing the Hessian of 00573 at the point specified by \a x 00574 (where \c n is the domain dimension for this ADFun<Base> object). 00575 */ 00576 template <class Base> 00577 template <class VectorBase, class VectorSet> 00578 VectorBase ADFun<Base>::SparseHessian( 00579 const VectorBase& x, const VectorBase& w, const VectorSet& p 00580 ) 00581 { size_t i, j, k; 00582 00583 size_t n = Domain(); 00584 VectorBase hes(n * n); 00585 00586 CPPAD_ASSERT_KNOWN( 00587 size_t(x.size()) == n, 00588 "SparseHessian: size of x not equal domain size for f." 00589 ); 00590 00591 typedef typename VectorSet::value_type Set_type; 00592 typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type; 00593 00594 // initialize the return value as zero 00595 Base zero(0); 00596 for(i = 0; i < n; i++) 00597 for(j = 0; j < n; j++) 00598 hes[i * n + j] = zero; 00599 00600 // arguments to SparseHessianCompute 00601 Pattern_type s; 00602 CppAD::vector<size_t> row; 00603 CppAD::vector<size_t> col; 00604 sparse_hessian_work work; 00605 bool transpose = false; 00606 sparsity_user2internal(s, p, n, n, transpose); 00607 k = 0; 00608 for(i = 0; i < n; i++) 00609 { s.begin(i); 00610 j = s.next_element(); 00611 while( j != s.end() ) 00612 { row.push_back(i); 00613 col.push_back(j); 00614 k++; 00615 j = s.next_element(); 00616 } 00617 } 00618 size_t K = k; 00619 VectorBase H(K); 00620 00621 // now we have folded this into the following case 00622 SparseHessianCompute(x, w, s, row, col, H, work); 00623 00624 // now set the non-zero return values 00625 for(k = 0; k < K; k++) 00626 hes[ row[k] * n + col[k] ] = H[k]; 00627 00628 return hes; 00629 } 00630 /*! 00631 Compute a sparse Hessian 00632 00633 The C++ source code coresponding to this operation is 00634 \verbatim 00635 hes = SparseHessian(x, w) 00636 \endverbatim 00637 00638 00639 \tparam Base 00640 is the base type for the recording that is stored in this 00641 ADFun<Base object. 00642 00643 \tparam VectorBase 00644 is a simple vector class with elements of the \a Base. 00645 00646 \param x [in] 00647 is a vector specifing the point at which to compute the Hessian. 00648 00649 \param w [in] 00650 The Hessian is computed for a weighted sum of the components 00651 of the function corresponding to this ADFun<Base> object. 00652 The argument \a w specifies the weights for each component. 00653 It must have size equal to the range dimension for this ADFun<Base> object. 00654 00655 \return 00656 Will be a vector of size \c n * n containing the Hessian of 00657 at the point specified by \a x 00658 (where \c n is the domain dimension for this ADFun<Base> object). 00659 */ 00660 template <class Base> 00661 template <class VectorBase> 00662 VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w) 00663 { size_t i, j, k; 00664 typedef CppAD::vectorBool VectorBool; 00665 00666 size_t m = Range(); 00667 size_t n = Domain(); 00668 00669 // determine the sparsity pattern p for Hessian of w^T F 00670 VectorBool r(n * n); 00671 for(j = 0; j < n; j++) 00672 { for(k = 0; k < n; k++) 00673 r[j * n + k] = false; 00674 r[j * n + j] = true; 00675 } 00676 ForSparseJac(n, r); 00677 // 00678 VectorBool s(m); 00679 for(i = 0; i < m; i++) 00680 s[i] = w[i] != 0; 00681 VectorBool p = RevSparseHes(n, s); 00682 00683 // compute sparse Hessian 00684 return SparseHessian(x, w, p); 00685 } 00686 00687 } // END_CPPAD_NAMESPACE 00688 # endif