Prev Next a_fun_hessian_xam.cpp Headings

@(@\newcommand{\B}[1]{ {\bf #1} } \newcommand{\R}[1]{ {\rm #1} }@)@
C++: Dense Hessian Using AD: Example and Test
# include <cstdio>
# include <string>
# include <cppad/swig/cppad_swig.hpp>

bool a_fun_hessian_xam(void) {
     using cppad_swig::a_double;
     using cppad_swig::vec_bool;
     using cppad_swig::vec_int;
     using cppad_swig::vec_double;
     using cppad_swig::vec_a_double;
     using cppad_swig::a_fun;
     using cppad_swig::sparse_rc;
     using cppad_swig::sparse_rcv;
     using cppad_swig::sparse_jac_work;
     using cppad_swig::sparse_hes_work;
     using std::string;
     //
     // initialize return variable
     bool ok = true;
     //------------------------------------------------------------------------
     // number of dependent and independent variables
     int n_dep = 1;
     int n_ind = 3;
     //
     // create the independent variables ax
     vec_double x = cppad_swig::vec_double(n_ind);
     for(int i = 0; i < n_ind ; i++) {
          x[i] = i + 2.0;
     }
     vec_a_double ax = cppad_swig::independent(x);
     //
     // create dependent variables ay with ay0 = ax_0 * ax_1 * ax_2
     a_double ax_0 = ax[0];
     a_double ax_1 = ax[1];
     a_double ax_2 = ax[2];
     vec_a_double ay = cppad_swig::vec_a_double(n_dep);
     ay[0] = ax_0 * ax_1 * ax_2;
     //
     // define af corresponding to f(x) = x_0 * x_1 * x_2
     a_fun af = cppad_swig::a_fun(ax, ay);
     //
     // g(x) = w_0 * f_0 (x) = f(x)
     vec_double w = cppad_swig::vec_double(n_dep);
     w[0] = 1.0;
     //
     // compute Hessian
     vec_double fpp = af.hessian(x, w);
     //
     //          [ 0.0 , x_2 , x_1 ]
     // f''(x) = [ x_2 , 0.0 , x_0 ]
     //          [ x_1 , x_0 , 0.0 ]
     ok = ok && fpp[0 * n_ind + 0] == 0.0 ;
     ok = ok && fpp[0 * n_ind + 1] == x[2] ;
     ok = ok && fpp[0 * n_ind + 2] == x[1] ;
     //
     ok = ok && fpp[1 * n_ind + 0] == x[2] ;
     ok = ok && fpp[1 * n_ind + 1] == 0.0 ;
     ok = ok && fpp[1 * n_ind + 2] == x[0] ;
     //
     ok = ok && fpp[2 * n_ind + 0] == x[1] ;
     ok = ok && fpp[2 * n_ind + 1] == x[0] ;
     ok = ok && fpp[2 * n_ind + 2] == 0.0 ;
     //
     return( ok );
}

Input File: build/lib/example/cplusplus/a_fun_hessian_xam.cpp