Prev Next atomic_four_lin_ode_base_solver.hpp

@(@\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}} }@)@This is cppad-20221105 documentation. Here is a link to its current documentation .
Atomic Multiply Base Matrices: Example Implementation

Syntax
lin_ode.base_solver(
    r, step, pattern, transpose, x, y
)

Prototype 

template <class Base>
void atomic_lin_ode<Base>::base_solver(
    const Base&                    r         ,
    const Base&                    step      ,
    const sparse_rc&               pattern   ,
    const bool&                    transpose ,
    const CppAD::vector<Base>&     x         ,
    CppAD::vector<Base>&           y         )

Notation
We use the notation: call_id r pattern transpose nnz , row , col , x , n , A(x) , b(x) , y(x) , m , vk(x)

Rosen34
This example uses one step of rosen34 ODE solver. Any initial value ODE solver, with any number of steps, could be used.

Source 

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

namespace CppAD { // BEGIN_CPPAD_NAMESPACE
//
// base_solver
// BEGIN_PROTOTYPE
template <class Base>
void atomic_lin_ode<Base>::base_solver(
    const Base&                    r         ,
    const Base&                    step      ,
    const sparse_rc&               pattern   ,
    const bool&                    transpose ,
    const CppAD::vector<Base>&     x         ,
    CppAD::vector<Base>&           y         )
// END_PROTOTYPE
{
    class Fun {
    private:
        const sparse_rc&           pattern_;
        const bool&                transpose_;
        const CppAD::vector<Base>& x_;
    public:
        Fun(
            const sparse_rc&           pattern   ,
            const bool&                transpose ,
            const CppAD::vector<Base>& x         )
        : pattern_(pattern), transpose_(transpose), x_(x)
        { }
        void Ode(
            const Base&                t ,
            const CppAD::vector<Base>& z ,
            CppAD::vector<Base>&       f )
        {   size_t m   = z.size();
            size_t nnz = pattern_.nnz();
            CPPAD_ASSERT_UNKNOWN( f.size() == m );
            CPPAD_ASSERT_UNKNOWN( x_.size() == nnz + m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nr() == m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nc() == m );
            //
            for(size_t i = 0; i < m; ++i)
                f[i] = Base(0);
            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);
                f[i] += x_[k] * z[j];
            }
        }
        void Ode_ind(
            const Base&                t   ,
            const CppAD::vector<Base>& z   ,
            CppAD::vector<Base>&       f_t )
        {   size_t m   = z.size();
# ifndef NDEBUG
            size_t nnz = pattern_.nnz();
            CPPAD_ASSERT_UNKNOWN( f_t.size() == m );
            CPPAD_ASSERT_UNKNOWN( x_.size() == nnz + m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nr() == m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nc() == m );
# endif
            //
            for(size_t i = 0; i < m; ++i)
                f_t[i] = Base(0);
        }
        void Ode_dep(
            const Base&                t   ,
            const CppAD::vector<Base>& z   ,
            CppAD::vector<Base>&       f_x )
        {   size_t m   = z.size();
            size_t nnz = pattern_.nnz();
            CPPAD_ASSERT_UNKNOWN( f_x.size() == m * m );
            CPPAD_ASSERT_UNKNOWN( x_.size() == nnz + m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nr() == m );
            CPPAD_ASSERT_UNKNOWN( pattern_.nc() == m );
            //
            for(size_t i = 0; i < m * m; ++i)
                f_x[i] = Base(0);
            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);
                f_x[i * m + j] = x_[k];
            }
        }
    };
    //
    // nnz
    size_t nnz = pattern.nnz();
    // m
    size_t m     = y.size();
    CPPAD_ASSERT_UNKNOWN( x.size() == nnz + m );
    //
    // fun
    Fun fun(pattern, transpose, x);
    //
    // y
    Base ti       = Base(0.0);
    Base tf       = r;
    size_t n_step = 1;
    if( step < abs(r) )
        n_step = size_t( Integer( abs(r) / step ) ) + 1;
    CppAD::vector<Base> zi(m), e(m);
    for(size_t j = 0; j < m; ++j)
        zi[j] = x[nnz + j];
    y = CppAD::Rosen34(fun, n_step, ti, tf, zi, e);
    return;
}
} // END_CPPAD_NAMESPACE
Input File: include/cppad/example/atomic_four/lin_ode/base_solver.hpp