\(\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}} }\)

adolc_sparse_jacobian.cpp

View page source

Adolc Speed: Sparse Jacobian

Specifications

See link_sparse_jacobian .

Implementation

// suppress conversion warnings before other includes
# include <cppad/wno_conversion.hpp>
//
# include <adolc/adolc.h>
# include <adolc/adolc_sparse.h>
# include <cppad/utility/vector.hpp>
# include <cppad/speed/uniform_01.hpp>
# include <cppad/speed/sparse_jac_fun.hpp>

// list of possible options
# include <map>
extern std::map<std::string, bool> global_option;

namespace {
   using CppAD::vector;
   typedef vector<size_t>    s_vector;
   typedef vector<double>    d_vector;
   typedef vector<adouble>   a_vector;
   void setup(
      // inputs
      int             tag     ,
      size_t          size    ,
      size_t          m       ,
      const s_vector& row     ,
      const s_vector& col     ,
      const d_vector& x       ,
      int*            options , // const but adolc want non-const arg
      // outupts
      size_t&         n_color ,
      int&            nnz     ,
      unsigned int*&  rind    ,
      unsigned int*&  cind    ,
      double*&        values  )
   {  // independent variables
      CPPAD_ASSERT_UNKNOWN( size = x.size() );
      int keep = 0; // keep forward mode results
      trace_on(tag, keep);
      size_t n = size;
      a_vector a_x(n);
      for(size_t j = 0; j < n; ++j)
         a_x[j] <<= x[j];
      //
      // dependent variables
      a_vector a_y(m);
      //
      // AD computation of f(x)
      size_t order = 0;
      CppAD::sparse_jac_fun<adouble>(m, n, a_x, row, col, order, a_y);
      //
      // create function object f : x -> y
      double yi;
      for(size_t i = 0; i < m; i++)
         a_y[i] >>= yi;
      trace_off();
      //
      // null pointers for recalculation of sparsity pattern
      free(rind);
      free(cind);
      free(values);
      rind   = nullptr;
      cind   = nullptr;
      values = nullptr;
      //
      // Retrieve n_color using undocumented feature of sparsedrivers.cpp
      int same_pattern = 0;
      n_color = sparse_jac(tag, int(m), int(n),
         same_pattern, x.data(), &nnz, &rind, &cind, &values, options
      );
   }

}

bool link_sparse_jacobian(
   const std::string&               job      ,
   size_t                           size     ,
   size_t                           repeat   ,
   size_t                           m        ,
   const CppAD::vector<size_t>&     row      ,
   const CppAD::vector<size_t>&     col      ,
          CppAD::vector<double>&     x_return ,
          CppAD::vector<double>&     jacobian ,
          size_t&                    n_color  )
{
   // --------------------------------------------------------------------
   // check global options
   // Allow colpack true even though it is not used below because it is
   // true durng the adolc correctness tests.
   const char* valid[] = { "onetape", "optimize", "colpack"};
   size_t n_valid = sizeof(valid) / sizeof(valid[0]);
   typedef std::map<std::string, bool>::iterator iterator;
   //
   for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
   {  if( itr->second )
      {  bool ok = false;
         for(size_t i = 0; i < n_valid; i++)
            ok |= itr->first == valid[i];
         if( ! ok )
            return false;
      }
   }
   // -----------------------------------------------------
   static size_t  static_size     = 0;
   static int     static_nnz      = 0;
   unsigned int*  static_rind     = nullptr;
   unsigned int*  static_cind     = nullptr;
   double*        static_values   = nullptr;
   // -----------------------------------------------------
   int tag  = 0;
   //
   // options that control sparse_jac
   int        options[4];
   if( global_option["boolsparsity"] )
      options[0] = 1;  // sparsity by propagation of bit pattern
   else
      options[0] = 0;  // sparsity pattern by index domains
   options[1] = 0;  // 0 = safe mode, 1 = tight mode
   options[2] = 0;  // 0 = autodetect, 1 = forward, 2 = reverse
   options[3] = 0;  // 0 = column compression, 1 = row compression
   //
   // independent variiables
   size_t n = size;
   d_vector x(n);
   //
   // default value for n_color
   n_color = 0;
   //
   bool onetape = global_option["onetape"];
   // -----------------------------------------------------
   if( job == "setup" )
   {  // get a value for x
      CppAD::uniform_01(n, x);
      //
      // record the tape and run coloring problem
      options[2] = -1;
      setup(tag, size, m, row, col, x, options,
         n_color, static_nnz, static_rind, static_cind, static_values
      );
      static_size = size;
      //
      return true;
   }
   if( job == "teardown" )
   {  free(static_rind);
      free(static_cind);
      free(static_values);
      static_rind   = nullptr;
      static_cind   = nullptr;
      static_values = nullptr;
      return true;
   }
   // -----------------------------------------------------
   CPPAD_ASSERT_UNKNOWN( job == "run" );
   //
   while (repeat--)
   {  // choose a value for x
      CppAD::uniform_01(n, x);
      //
      if ( ! onetape )
      {  // retape and calculate jacobian
         options[2] = -1; // stop at sparsity pattern, return n_color
         setup(tag, size, m, row, col, x, options,
            n_color, static_nnz, static_rind, static_cind, static_values
         );
         options[2] = 0;
      }
      else
      {  if( size != static_size )
            CPPAD_ASSERT_UNKNOWN( size == static_size );
      }
      // calculate the jacobian at this x
      int same_pattern = 1;
      sparse_jac(
         tag, int(m), int(n), same_pattern, x.data(),
         &static_nnz, &static_rind, &static_cind, &static_values,
         options
      );
   }
   // --------------------------------------------------------------------
   // jacobian
   CPPAD_ASSERT_UNKNOWN( size_t(static_nnz) == row.size() );
   for(int ell = 0; ell < static_nnz; ell++)
   {  CPPAD_ASSERT_UNKNOWN( row[ell] == size_t(static_rind[ell]) );
      CPPAD_ASSERT_UNKNOWN( col[ell] == size_t(static_cind[ell]) );
      jacobian[ell] = static_values[ell];
   }
   // x_return
   for(size_t j = 0; j < n; j++)
      x_return[j] = x[j];
   //
   return true;
}