Prev

Next

Index-> contents reference index search external

Up-> CppAD speed speed_adolc adolc_ode.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_ode.cpp

Headings-> Specifications Implementation

// suppress conversion warnings before other includes
# include <cppad/wno_conversion.hpp>
//
# include <adolc/adolc.h>

# include <cppad/utility/vector.hpp>
# include <cppad/speed/ode_evaluate.hpp>
# include <cppad/speed/uniform_01.hpp>

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

bool link_ode(
     size_t                     size       ,
     size_t                     repeat     ,
     CppAD::vector<double>      &x         ,
     CppAD::vector<double>      &jac
)
{
     // speed test global option values
     if( global_option["atomic"] )
          return false;
     if( global_option["memory"] || global_option["optimize"] )
          return false;
     // -------------------------------------------------------------
     // setup
     assert( x.size() == size );
     assert( jac.size() == size * size );

     typedef CppAD::vector<adouble> ADVector;
     typedef CppAD::vector<double>  DblVector;

     size_t i, j;
     int tag    = 0;       // tape identifier
     int keep   = 0;       // do not keep forward mode results
     size_t p   = 0;       // use ode to calculate function values
     size_t n   = size;    // number of independent variables
     size_t m   = n;       // number of dependent variables
     ADVector  X(n), Y(m); // independent and dependent variables
     DblVector f(m);       // function value

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

     // raw memory for use with adolc
     size_min = n;
     double *x_raw   = thread_alloc::create_array<double>(size_min, size_out);
     size_min = m * n;
     double *jac_raw = thread_alloc::create_array<double>(size_min, size_out);
     size_min = m;
     double **jac_ptr = thread_alloc::create_array<double*>(size_min, size_out);
     for(i = 0; i < m; i++)
          jac_ptr[i] = jac_raw + i * n;

     // -------------------------------------------------------------
     if( ! global_option["onetape"] ) while(repeat--)
     {     // choose next x value
          uniform_01(n, x);

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

          // evaluate function
          CppAD::ode_evaluate(X, p, Y);

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

          // evaluate the Jacobian
          for(j = 0; j < n; j++)
               x_raw[j] = x[j];
          jacobian(tag, m, n, x_raw, jac_ptr);
     }
     else
     {     // choose next x value
          uniform_01(n, x);

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

          // evaluate function
          CppAD::ode_evaluate(X, p, Y);

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

          while(repeat--)
          {     // get next argument value
               uniform_01(n, x);
               for(j = 0; j < n; j++)
                    x_raw[j] = x[j];

               // evaluate jacobian
               jacobian(tag, m, n, x_raw, jac_ptr);
          }
     }
     // convert return value to a simple vector
     for(i = 0; i < m; i++)
     {     for(j = 0; j < n; j++)
               jac[i * n + j] = jac_ptr[i][j];
     }
     // ----------------------------------------------------------------------
     // tear down
     thread_alloc::delete_array(x_raw);
     thread_alloc::delete_array(jac_raw);
     thread_alloc::delete_array(jac_ptr);

     return true;
}