Prev Next adolc_sparse_hessian.cpp

Adolc Speed: Sparse Hessian

Operation Sequence
Note that the operation sequence depends on the vectors i and j. Hence we use a different ADFun object for each choice of i and j.

link_sparse_hessian
Routine that computes the gradient of determinant using Adolc:
# include <cppad/vector.hpp>
# include <cppad/speed/uniform_01.hpp>
# include <cppad/track_new_del.hpp>
# include <cppad/speed/sparse_evaluate.hpp>

# include <adolc/adouble.h>
# include <adolc/drivers/drivers.h>

bool link_sparse_hessian(
     size_t                     repeat   , 
     CppAD::vector<double>     &x_arg    ,
     CppAD::vector<size_t>     &i        ,
     CppAD::vector<size_t>     &j        ,
     CppAD::vector<double>     &h        )
{
     // -----------------------------------------------------
     // setup
     size_t k, m;
     size_t order = 0;         // derivative order corresponding to function
     size_t tag  = 0;          // tape identifier
     size_t keep = 1;          // keep forward mode results in buffer
     size_t n = x_arg.size();  // number of independent variables
     size_t ell = i.size();    // number of indices in i and j
     double f;                 // function value

     typedef CppAD::vector<double>  DblVector;
     typedef CppAD::vector<adouble> ADVector;
     typedef CppAD::vector<size_t>  SizeVector;

     ADVector   X(n);    // AD domain space vector
     double       *x;    // double domain space vector
     double      **H;    // Hessian 
     ADVector   Y(1);    // AD range space value
     DblVector tmp(2 * ell);       // double temporary vector

     x = 0;
     x = CPPAD_TRACK_NEW_VEC(n, x);
     H = 0;
     H = CPPAD_TRACK_NEW_VEC(n, H);
     for(k = 0; k < n; k++)
     {    H[k] = 0;
          H[k] = CPPAD_TRACK_NEW_VEC(n, H[k]);
     }

     // choose a value for x 
     CppAD::uniform_01(n, x);
     for(k = 0; k < n; k++)
          x_arg[k] = x[k];

     // ------------------------------------------------------
     while(repeat--)
     {
          // get the next set of indices
          CppAD::uniform_01(2 * ell, tmp);
          for(k = 0; k < ell; k++)
          {    i[k] = size_t( n * tmp[k] );
               i[k] = std::min(n-1, i[k]);
               //
               j[k] = size_t( n * tmp[k + ell] );
               j[k] = std::min(n-1, j[k]);
          }

          // declare independent variables
          trace_on(tag, keep);
          for(k = 0; k < n; k++)
               X[k] <<= x[k];

          // AD computation of f(x)
          CppAD::sparse_evaluate(X, i, j, order, Y);

          // create function object f : X -> Y
          Y[0] >>= f;
          trace_off();

          // evaluate and return the hessian of f
          hessian(int(tag), int(n), x, H);
     }
     for(k = 0; k < n; k++)
     {    for(m = 0; m <= k; m++)
          {    h[ k * n + m] = H[k][m];
               h[ m * n + k] = H[k][m];
          }
          CPPAD_TRACK_DEL_VEC(H[k]);
     }
     CPPAD_TRACK_DEL_VEC(H);
     CPPAD_TRACK_DEL_VEC(x);
     return true;
}

Input File: speed/adolc/sparse_hessian.cpp