Atomic Linear ODE Hessian Sparsity Pattern: Example Implementation

atomic_four_lin_ode_hes_sparsity.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 Linear ODE Hessian Sparsity Pattern: Example Implementation
Purpose
The hes_sparsity routine overrides the virtual functions used by the atomic_four base class for Hessian sparsity calculations; see hes_sparsity .

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

wk(x)
Because we are using the Rosen34 solver, our actual sequence of operations is only fourth order accurate. So it suffices to compute the sparsity pattern for @[@ \tilde{y} (x) = \sum_{k=0}^4 v^k (x) @]@ Note that the factor @(@ r / k @)@, in the definition of @(@ v^k (x) @)@, is constant (with respect to the variables). Hence it suffices to compute the sparsity pattern for @[@ h (x) = \sum_{k=0}^4 w^k (x) @]@ where @(@ w^0 (x) = b(x) @)@ and for @(@ k = 1, 2, \ldots @)@, @(@ w^k (x) = A(x) w^{k-1} (x) @)@.

Example
The file atomic_four_lin_ode_sparsity.cpp contains an example and test using this operator.

Source


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

namespace CppAD { // BEGIN_CPPAD_NAMESPACE
//
// hes_sparsity override
template <class Base>
bool atomic_lin_ode<Base>::hes_sparsity(
    size_t                                         call_id      ,
    const CppAD::vector<bool>&                     ident_zero_x ,
    const CppAD::vector<bool>&                     select_x     ,
    const CppAD::vector<bool>&                     select_y     ,
    CppAD::sparse_rc< CppAD::vector<size_t> >&     pattern_out  )
{   //
    // adouble
    typedef AD<double> adouble;
    //
    // pattern_A, transpose, nnz
    Base      r;
    Base      step;
    sparse_rc pattern_A;
    bool      transpose;
    get(call_id, r, step, pattern_A, transpose);
    size_t nnz = pattern_A.nnz();
    //
    // m, n
    size_t m = select_y.size();
    size_t n = select_x.size();
    //
    // au
    vector<adouble> au(n);
    for(size_t j = 0; j < n; ++j)
        au[j] = 1.0;
    Independent( au );
    //
    // ax
    vector<adouble> ax(n);
    for(size_t j = 0; j < n; ++j)
        if( ident_zero_x[j] )
            ax[j] = 0.0;
        else
            ax[j] = au[j];
    //
    // aw
    // aw = w^0 (x)
    vector<adouble> aw(m);
    for(size_t i = 0; i < m; ++i)
        aw[i] = ax[nnz + i];
    //
    // ah = w^0 (x)
    vector<adouble> ah = aw;
    //
    // ah = sum_k=0^4 w^k (x)
    vector<adouble> awk(m);
    for(size_t k = 1; k < 5; ++k)
    {   // aw = w^{k-1} (x)
        //
        // awk = w^k (x)
        for(size_t i = 0; i < m; ++i)
            awk[i] = 0.0;
        for(size_t p = 0; p < nnz; ++p)
        {   //
            // i, j
            size_t i = pattern_A.row()[p];
            size_t j = pattern_A.col()[p];
            if( transpose )
                std::swap(i, j);
            //
            // awk
            awk[i] += ax[p] * aw[j];
        }
        //
        // ah = ah + w^k(x)
        for(size_t i = 0; i < m; ++i)
            ah[i] += awk[i];
        //
        // aw = w^k (x)
        aw = awk;
    }
    //
    // h
    ADFun<double> h;
    h.Dependent(au, ah);
    //
    // pattern_out
    // can use h.for_hes_sparsity or h.rev_hes_sparsity
    bool internal_bool = false;
    h.for_hes_sparsity(select_x, select_y, internal_bool, pattern_out);
    //
    return true;
}
} // END_CPPAD_NAMESPACE

Input File: include/cppad/example/atomic_four/lin_ode/hes_sparsity.hpp