Prev Next sparse_hes_pattern_xam.cpp Headings

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

bool sparse_hes_pattern_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 = 3;
     //
     // create the independent variables ax
     vec_double x = cppad_swig::vec_double(n);
     for(int i = 0; i < n ; i++) {
          x[i] = i + 2.0;
     }
     vec_a_double ax = cppad_swig::independent(x);
     //
     // create dependent variables ay with ay[i] = ax[j] * ax[i]
     // where i = mod(j + 1, n)
     vec_a_double ay = cppad_swig::vec_a_double(n);
     for(int j = 0; j < n ; j++) {
          int i = j+1;
          if( i >= n  ) {
               i = i - n;
          }
          a_double ay_i = ax[i] * ax[j];
          ay[i] = ay_i;
     }
     //
     // define af corresponding to f(x)
     a_fun af = cppad_swig::a_fun(ax, ay);
     //
     // Set select_d (domain) to all true, initial select_r (range) to all false
     vec_bool select_d = cppad_swig::vec_bool(n);
     vec_bool select_r = cppad_swig::vec_bool(n);
     for(int i = 0; i < n; i++) {
          select_d[i] = true;
          select_r[i] = false;
     }
     //
     // only select component 0 of the range function
     // f_0 (x) = x_0 * x_{n-1}
     select_r[0] = true;
     //
     // loop over forward and reverse mode
     for(int mode = 0; mode < 2; mode++) {
          sparse_rc pat_out = cppad_swig::sparse_rc();
          if( mode == 0  ) {
               af.for_hes_sparsity(select_d, select_r, pat_out);
          }
          if( mode == 1  ) {
               af.rev_hes_sparsity(select_d, select_r, pat_out);
          }
          //
          // check that result is sparsity pattern for Hessian of f_0 (x)
          ok = ok && pat_out.nnz() == 2 ;
          vec_int row = pat_out.row();
          vec_int col = pat_out.col();
          for(int k = 0; k < 2; k++) {
               int r = row[k];
               int c = col[k];
               if( r <= c  ) {
                    ok = ok && r == 0;
                    ok = ok && c == n-1;
               }
               if( r >= c  ) {
                    ok = ok && r == n-1;
                    ok = ok && c == 0;
               }
          }
     }
     //
     return( ok );
}

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