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

atomic_four_lin_ode_forward.hpp

View page source

Atomic Linear ODE Forward Mode: Example Implementation

Purpose

The forward routine overrides the virtual functions used by the atomic_four base; see forward .

Theory

Suppose we are given Taylor coefficients \(x^0\), \(x^1\), for x . The zero order Taylor coefficient for z(t, x) solves the following initial value ODE:

\[z_t^0 (t) = A^0 z(t) \W{,} z^0 (0) = b^0\]

Note that \(A^0\) and \(b^0\) are just certain components of \(x^0\); see A(x) and b(x) . The first order Taylor coefficient for \(z(t, x)\) solves the following initial value ODE:

\[z_t^1 (t) = A^0 z^1 (t) + A^1 z^0 (t) \W{,} z^1 (0) = b^1\]

Note that \(A^1\) and \(c^1\) are just certain components of \(x^1\). We can solve for \(z^1 (t)\) using the following extended initial value ODE:

\[\begin{split}\left[ \begin{array}{c} z^0_t (t, x) \\ z^1_t (t, x) \end{array} \right] = \left[ \begin{array}{cc} A^0 & 0 \\ A^1 & A^0 \end{array} \right] \left[ \begin{array}{c} z^0 (t, x) \\ z^1 (t, x) \end{array} \right] \; , \; \left[ \begin{array}{c} z^0 (0, x) \\ z^1 (0, x) \end{array} \right] = \left[ \begin{array}{c} b^0 \\ b^1 \end{array} \right]\end{split}\]

extend_ode

The extended system above is created form the original system by the extend_ode function defined below:

Source

# include <cppad/example/atomic_four/lin_ode/lin_ode.hpp>

namespace CppAD { // BEGIN_CPPAD_NAMESPACE
// ----------------------------------------------------------------------------
// forward override for Base atomic linear ODE
template <class Base>
bool atomic_lin_ode<Base>::forward(
   size_t                                     call_id     ,
   const CppAD::vector<bool>&                 select_y    ,
   size_t                                     order_low   ,
   size_t                                     order_up    ,
   const CppAD::vector<Base>&                 taylor_x    ,
   CppAD::vector<Base>&                       taylor_y    )
{
   // order_up
   if( order_up > 1 )
      return false;
   //
   // r, pattern, transpose, nnz
   Base      r;
   Base      step;
   sparse_rc pattern;
   bool      transpose;
   get(call_id, r, step, pattern, transpose);
# ifndef NDEBUG
   size_t nnz = pattern.nnz();
# endif
   //
   // q
   size_t q = order_up + 1;
   //
   // m
   CPPAD_ASSERT_UNKNOWN( taylor_y.size() % q == 0 );
   size_t m = taylor_y.size() / q;
   CPPAD_ASSERT_UNKNOWN( pattern.nr() == m );
   CPPAD_ASSERT_UNKNOWN( pattern.nc() == m );
   //
   // taylor_x
   CPPAD_ASSERT_UNKNOWN( taylor_x.size() == (nnz + m) * q );
   //
   // taylor_y
   if( order_up == 0 )
      base_solver(r, step, pattern, transpose, taylor_x, taylor_y);
   else
   {  CPPAD_ASSERT_UNKNOWN( order_up == 1 );
      //
      // pattern_extend, x_extend
      sparse_rc           pattern_extend;
      CppAD::vector<Base> x_extend;
      extend_ode(pattern, transpose, taylor_x, pattern_extend, x_extend);
      //
      // y_extend
      size_t m_extend = pattern_extend.nr();
      CppAD::vector<Base> y_extend(m_extend);
      bool transpose_extend = false;
      base_solver(
         r, step, pattern_extend, transpose_extend, x_extend, y_extend
      );
      //
      // taylor_y
      if( order_low == 0 )
      {  for(size_t i = 0; i < m; ++i)
            taylor_y[i * q + 0] = y_extend[i];
      }
      for(size_t i = 0; i < m; ++i)
         taylor_y[i * q + 1] = y_extend[m + i];
   }
   //
   return true;
}
// ----------------------------------------------------------------------------
// forward override for AD<Base> atomic linear ODE
template <class Base>
bool atomic_lin_ode<Base>::forward(
   size_t                                     call_id     ,
   const CppAD::vector<bool>&                 select_y    ,
   size_t                                     order_low   ,
   size_t                                     order_up    ,
   const CppAD::vector< CppAD::AD<Base> >&    ataylor_x   ,
   CppAD::vector< CppAD::AD<Base> >&          ataylor_y   )
{  //
   // ad_base
   typedef CppAD::AD<Base> ad_base;
   //
   // order_up
   if( order_up > 1 )
      return false;
   //
   // r, pattern, transpose, nnz
   Base            r;
   Base            step;
   sparse_rc       pattern;
   bool            transpose;
   get(call_id, r, step, pattern, transpose);
# ifndef NDEBUG
   size_t nnz = pattern.nnz();
# endif
   //
   // q
   size_t q = order_up + 1;
   //
   // m
   CPPAD_ASSERT_UNKNOWN( ataylor_y.size() % q == 0 );
   size_t m = ataylor_y.size() / q;
   //
   // ataylor_x
   CPPAD_ASSERT_UNKNOWN( ataylor_x.size() == (nnz + m) * q );
   //
   // ataylor_y
   if( order_up == 0 )
      (*this)(call_id, ataylor_x, ataylor_y);
   else
   {  CPPAD_ASSERT_UNKNOWN( order_up == 1 );
      //
      // pattern_extend, x_extend
      sparse_rc              pattern_extend;
      CppAD::vector<ad_base> ax_extend;
      extend_ode(pattern, transpose, ataylor_x, pattern_extend, ax_extend);
      //
      // call_id_2
      bool transpose_extend = false;
      size_t call_id_2 = set(r, step, pattern_extend, transpose_extend);
      //
      // ay_extend
      size_t m_extend = pattern_extend.nr();
      CppAD::vector<ad_base> ay_extend(m_extend);
      (*this)(call_id_2, ax_extend, ay_extend);
      //
      // ataylor_y
      if( order_low == 0 )
      {  for(size_t i = 0; i < m; ++i)
            ataylor_y[i * q + 0] = ay_extend[i];
      }
      for(size_t i = 0; i < m; ++i)
         ataylor_y[i * q + 1] = ay_extend[m + i];
   }
   //
   return true;
}
// ---------------------------------------------------------------------------
// extend_ode
/*
[   z_t^0 (t, x^0 )    ] = [ A^0   0    ] * [   z^0 (t, x^0 )    ]
[ z_t^1 (t, x^0, x^1 ) ]   [ A^1   A^0  ]   [ z^1 (t, x^0, x^1 ) ]

[   z^0 (0, x^0 )    ] = [ b^0 ]
[ z^1 (0, x^0, x^1 ) ]   [ b^1 ]
*/
template <class Base>
template <class Float>
void atomic_lin_ode<Base>::extend_ode(
   const CppAD::sparse_rc< CppAD::vector<size_t> >&    pattern         ,
   const bool&                                         transpose       ,
   const CppAD::vector<Float>&                         taylor_x        ,
   CppAD::sparse_rc< CppAD::vector<size_t> >&          pattern_extend  ,
   CppAD::vector<Float>&                               x_extend        )
{  //
   // m
   size_t m   = pattern.nr();
   size_t nnz = pattern.nnz();
   CPPAD_ASSERT_UNKNOWN( pattern.nc() == m );
   //
   // q
   size_t q = 2;
   CPPAD_ASSERT_UNKNOWN( taylor_x.size() == q * (nnz + m) );
   //
   // m_extend
   size_t m_extend   = 2 * m;
   size_t nnz_extend = 3 * nnz;
   //
   // pattern_extend, x_extend
   pattern_extend.resize(m_extend, m_extend, nnz_extend);
   x_extend.resize(nnz_extend + m_extend);
   for(size_t k = 0; k < nnz; ++k)
   {  size_t i = pattern.row()[k];
      size_t j = pattern.col()[k];
      if( transpose )
         std::swap(i, j);
      //
      // A^0_ij
      Float A0ij = taylor_x[k * q + 0];
      //
      // A^1_ij
      Float A1ij = taylor_x[k * q + 1];
      //
      // upper diagonal
      pattern_extend.set(3 * k + 0, i, j);
      x_extend[3 * k + 0] = A0ij;
      //
      // lower left
      pattern_extend.set(3 * k + 1, m + i, j);
      x_extend[3 * k + 1] = A1ij;
      //
      // lower diagonal
      pattern_extend.set(3 * k + 2, m + i, m + j);
      x_extend[3 * k + 2] = A0ij;
   }
   for(size_t i = 0; i < m; ++i)
   {  // b^0_i
      x_extend[nnz_extend + i]     = taylor_x[ (nnz + i) * q + 0 ];
      // b^1_i
      x_extend[nnz_extend + m + i] = taylor_x[ (nnz + i) * q + 1 ];
   }
   return;
}
} // END_CPPAD_NAMESPACE