Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse sparse_jac sparse_jac_xam.cpp

cppad_swig-> testing swig_xam.i library whats_new_2018

library-> module a_double vector a_fun sparse error xam.m4

sparse-> sparse_rc sparse_rcv jac_sparsity hes_sparsity sparse_jac sparse_hes

sparse_jac-> sparse_jac_xam.cpp sparse_jac_xam.m sparse_jac_xam.pm sparse_jac_xam.py

sparse_jac_xam.cpp

Headings

# include <cstdio>
# include <string>
# include <cppad/swig/cppad_swig.hpp>

bool sparse_jac_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;
     // one
     a_double aone = cppad_swig::a_double(1.0);
     //
     // 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] = (j+1) * ax[j]
     // 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 aj = cppad_swig::a_double(j);
          a_double ay_i = (aj + aone) * ax[j];
          ay[i] = ay_i;
     }
     //
     // define af corresponding to f(x)
     a_fun af = cppad_swig::a_fun(ax, ay);
     //
     // sparsity pattern for identity matrix
     sparse_rc pat_eye = cppad_swig::sparse_rc();
     pat_eye.resize(n, n, n);
     for(int k = 0; k < n; k++) {
          pat_eye.put(k, k, k);
     }
     //
     // sparsity pattern for the Jacobian
     sparse_rc pat_jac = cppad_swig::sparse_rc();
     af.for_jac_sparsity(pat_eye, pat_jac);
     //
     // loop over forward and reverse mode
     for(int mode = 0; mode < 2; mode++) {
          // compute all possibly non-zero entries in Jacobian
          sparse_rcv subset = cppad_swig::sparse_rcv(pat_jac);
          // work space used to save time for multiple calls
          sparse_jac_work work = cppad_swig::sparse_jac_work();
          if( mode == 0  ) {
               af.sparse_jac_for(subset, x, pat_jac, work);
          }
          if( mode == 1  ) {
               af.sparse_jac_rev(subset, x, pat_jac, work);
          }
          //
          // check result
          ok = ok && n == subset.nnz();
          vec_int col_major = subset.col_major();
          vec_int row = subset.row();
          vec_int col = subset.col();
          vec_double val = subset.val();
          for(int k = 0; k < n; k++) {
               int ell = col_major[k];
               int r = row[ell];
               int c = col[ell];
               double v = val[ell];
               int i = c+1;
               if( i >=  n  ) {
                    i = i - n;
               }
               ok = ok && c == k;
               ok = ok && r == i;
               ok = ok && v == c + 1.0;
          }
     }
     //
     return( ok );
}