Prev

Next

Index-> contents reference index search external

Up-> CppAD speed speed_adolc adolc_sparse_jacobian.cpp

CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix

speed-> speed_main speed_utility speed_double speed_adolc speed_cppad speed_fadbad speed_sacado

speed_adolc-> adolc_det_minor.cpp adolc_det_lu.cpp adolc_mat_mul.cpp adolc_ode.cpp adolc_poly.cpp adolc_sparse_hessian.cpp adolc_sparse_jacobian.cpp adolc_alloc_mat

adolc_sparse_jacobian.cpp

Headings-> Specifications 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;

bool link_sparse_jacobian(
     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_sweep  )
{
     if( global_option["atomic"] || (! global_option["colpack"]) )
          return false;
     if( global_option["memory"] || global_option["optimize"] )
          return false;
     // -----------------------------------------------------
     // setup
     typedef unsigned int*    SizeVector;
     typedef double*          DblVector;
     typedef adouble          ADScalar;
     typedef ADScalar*        ADVector;

     size_t i, j, k;            // temporary indices
     size_t n = size;           // number of independent variables
     size_t order = 0;          // derivative order corresponding to function

     // set up for thread_alloc memory allocator (fast and checks for leaks)
     using CppAD::thread_alloc; // the allocator
     size_t capacity;           // capacity of an allocation

     // tape identifier
     int tag  = 0;
     // AD domain space vector
     ADVector a_x = thread_alloc::create_array<ADScalar>(n, capacity);
     // AD range space vector
     ADVector a_y = thread_alloc::create_array<ADScalar>(m, capacity);
     // argument value in double
     DblVector x = thread_alloc::create_array<double>(n, capacity);
     // function value in double
     DblVector y = thread_alloc::create_array<double>(m, capacity);


     // 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; // see changing to -1 and back to 0 below
     options[3] = 0; // (0 = column compression, 1 = row compression)

     // structure that holds some of the work done by sparse_jac
     int        nnz;                   // number of non-zero values
     SizeVector rind   = CPPAD_NULL;   // row indices
     SizeVector cind   = CPPAD_NULL;   // column indices
     DblVector  values = CPPAD_NULL;   // Jacobian values

     // choose a value for x
     CppAD::uniform_01(n, x);

     // declare independent variables
     int keep = 0; // keep forward mode results
     trace_on(tag, keep);
     for(j = 0; j < n; j++)
          a_x[j] <<= x[j];

     // AD computation of f (x)
     CppAD::sparse_jac_fun<ADScalar>(m, n, a_x, row, col, order, a_y);

     // create function object f : x -> y
     for(i = 0; i < m; i++)
          a_y[i] >>= y[i];
     trace_off();

     // Retrieve n_sweep using undocumented feature of sparsedrivers.cpp
     int same_pattern = 0;
     options[2]       = -1;
     n_sweep = sparse_jac(tag, int(m), int(n),
          same_pattern, x, &nnz, &rind, &cind, &values, options
     );
     options[2]       = 0;
     // ----------------------------------------------------------------------
     if( ! global_option["onetape"] ) while(repeat--)
     {     // choose a value for x
          CppAD::uniform_01(n, x);

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

          // AD computation of f (x)
          CppAD::sparse_jac_fun<ADScalar>(m, n, a_x, row, col, order, a_y);

          // create function object f : x -> y
          for(i = 0; i < m; i++)
               a_y[i] >>= y[i];
          trace_off();

          // is this a repeat call with the same sparsity pattern
          same_pattern = 0;

          // calculate the jacobian at this x
          rind   = CPPAD_NULL;
          cind   = CPPAD_NULL;
          values = CPPAD_NULL;
          sparse_jac(tag, int(m), int(n),
               same_pattern, x, &nnz, &rind, &cind, &values, options
          );
          // only needed last time through loop
          if( repeat == 0 )
          {     size_t K = row.size();
               for(int ell = 0; ell < nnz; ell++)
               {     i = size_t(rind[ell]);
                    j = size_t(cind[ell]);
                    for(k = 0; k < K; k++)
                    {     if( row[k]==i && col[k]==j )
                              jacobian[k] = values[ell];
                    }
               }
          }

          // free raw memory allocated by sparse_jac
          free(rind);
          free(cind);
          free(values);
     }
     else
     {     while(repeat--)
          {     // choose a value for x
               CppAD::uniform_01(n, x);

               // calculate the jacobian at this x
               sparse_jac(tag, int(m), int(n),
                    same_pattern, x, &nnz, &rind, &cind, &values, options
               );
               same_pattern = 1;
          }
          size_t K = row.size();
          for(int ell = 0; ell < nnz; ell++)
          {     i = size_t(rind[ell]);
               j = size_t(cind[ell]);
               for(k = 0; k < K; k++)
               {     if( row[k]==i && col[k]==j )
                         jacobian[k] = values[ell];
               }
          }

          // free raw memory allocated by sparse_jac
          free(rind);
          free(cind);
          free(values);
     }
     // --------------------------------------------------------------------
     // return argument
     for(j = 0; j < n; j++)
          x_return[j] = x[j];

     // tear down
     thread_alloc::delete_array(a_x);
     thread_alloc::delete_array(a_y);
     thread_alloc::delete_array(x);
     thread_alloc::delete_array(y);
     return true;
}