Prev Next sparse_hes_fun.hpp Headings

@(@\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}} }@)@
Source: sparse_hes_fun
# 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

Input File: omh/sparse_hes_fun.omh