CppAD: A C++ Algorithmic Differentiation Package  20130918
sparse_hessian.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines