$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
# ifndef CPPAD_SPARSE_HES_FUN_HPP # define CPPAD_SPARSE_HES_FUN_HPP # include <cppad/core/cppad_assert.hpp> # include <cppad/utility/check_numeric_type.hpp> # include <cppad/utility/vector.hpp> // following needed by gcc under fedora 17 so that exp(double) is defined # include <cppad/base_require.hpp> namespace CppAD { template <class Float, class FloatVector> void sparse_hes_fun( size_t n , const FloatVector& x , const CppAD::vector<size_t>& row , const CppAD::vector<size_t>& col , size_t p , FloatVector& fp ) { // check numeric type specifications CheckNumericType<Float>(); // check value of p CPPAD_ASSERT_KNOWN( p == 0 || p == 2, "sparse_hes_fun: p != 0 and p != 2" ); size_t K = row.size(); size_t i, j, k; if( p == 0 ) fp[0] = Float(0); else { for(k = 0; k < K; k++) fp[k] = Float(0); } // determine which diagonal entries are present in row[k], col[k] CppAD::vector<size_t> diagonal(n); for(i = 0; i < n; i++) diagonal[i] = K; // no diagonal entry for this row for(k = 0; k < K; k++) { if( row[k] == col[k] ) { CPPAD_ASSERT_UNKNOWN( diagonal[row[k]] == K ); // index of the diagonal entry diagonal[ row[k] ] = k; } } // determine which entries must be multiplied by a factor of two CppAD::vector<Float> factor(K); for(k = 0; k < K; k++) { factor[k] = Float(1); for(size_t k1 = 0; k1 < K; k1++) { bool reflected = true; reflected &= k != k1; reflected &= row[k] != col[k]; reflected &= row[k] == col[k1]; reflected &= col[k] == row[k1]; if( reflected ) factor[k] = Float(2); } } Float t; for(k = 0; k < K; k++) { i = row[k]; j = col[k]; t = exp( x[i] * x[j] ); switch(p) { case 0: fp[0] += t; break; case 2: if( i == j ) { // dt_dxi = 2.0 * xi * t fp[k] += ( Float(2) + Float(4) * x[i] * x[i] ) * t; } else { // dt_dxi = xj * t fp[k] += factor[k] * ( Float(1) + x[i] * x[j] ) * t; if( diagonal[i] != K ) { size_t ki = diagonal[i]; fp[ki] += x[j] * x[j] * t; } if( diagonal[j] != K ) { size_t kj = diagonal[j]; fp[kj] += x[i] * x[i] * t; } } break; } } } } # endif