Use of odeint with boost multi_array

[jannisehrlich]

Hi I am quite new to C++ and, thus, boost. In the end I need to integrate a 10 dimensional array with odeint, which is a private member of a class. I am testing the compatibility with a small example of the harmonic oscillator - which I attached (github somehow does not enable uploading a cpp file), compiling with g++ -Wall -fexceptions -std=c++11 -g -c -o Test_odeint.x Test_odeint.cpp When I used boost::multi_array<double, 1> v_array; it worked, but switching to boost::multi_array<double, 2> v_array according to the example code produces an error in the for_each subroutine, esp. for_each_unrolled - as far as I understand the error messages.

I encountered a similar problem directly using the multi_array class (commented in the attached file) where I received similar errors.

Test_odeint.cpp.txt

I could not find out, what I did wrong by reading the rather lengthy error messages:

/usr/local/include/boost/fusion/algorithm/iteration/detail/for_each.hpp|135| [ skipping 2 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]|

/home/jehrlich/programming/test_ode_solver/Test_ode_solver/main.cpp|316|required from here| /usr/local/include/boost/numeric/odeint/algebra/detail/for_each.hpp|47|error: no match for call to ‘(boost::numeric::odeint::detail::generic_rk_scale_sum<1ul, boost::numeric::odeint::default_operations, double, double>) (boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference)’

/usr/local/include/boost/fusion/algorithm/iteration/detail/for_each.hpp|135| [ skipping 2 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]

/home/jehrlich/programming/test_ode_solver/Test_ode_solver/main.cpp|316|required from here| /usr/local/include/boost/numeric/odeint/algebra/detail/for_each.hpp|55|error: no match for call to ‘(boost::numeric::odeint::detail::generic_rk_scale_sum<2ul, boost::numeric::odeint::default_operations, double, double>) (boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference)’

/usr/local/include/boost/fusion/algorithm/iteration/detail/for_each.hpp|135| [ skipping 2 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]

/home/jehrlich/programming/test_ode_solver/Test_ode_solver/main.cpp|316|required from here| /usr/local/include/boost/numeric/odeint/algebra/detail/for_each.hpp|64|error: no match for call to ‘(boost::numeric::odeint::detail::generic_rk_scale_sum<3ul, boost::numeric::odeint::default_operations, double, double>) (boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference)’

/usr/local/include/boost/fusion/algorithm/iteration/detail/for_each.hpp|135| [ skipping 2 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]

/home/jehrlich/programming/test_ode_solver/Test_ode_solver/main.cpp|316|required from here| /usr/local/include/boost/numeric/odeint/algebra/detail/for_each.hpp|73|error: no match for call to ‘(boost::numeric::odeint::detail::generic_rk_scale_sum<4ul, boost::numeric::odeint::default_operations, double, double>) (boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, const double*, mpl::size_t<2ul>, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::const_sub_array<double, 1ul, const double>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference, boost::iterators::detail::iterator_facade_base<boost::detail::multi_array::arrayiterator<double, double*, mpl::size_t<2ul>, boost::detail::multi_array::sub_array<double, 1ul>, boost::iterators::random_access_traversal_tag>, boost::multi_array<double, 1ul, std::allocator >, boost::iterators::random_access_traversal_tag, boost::detail::multi_array::sub_array<double, 1ul>, long int, false, false>::reference)’

Can you please give me a hint what I have to do to get it working?

[mariomulansky]

Quickly skimming through your code it seem your harm_osc class is missing a `const in the operator():

void operator() ( Statetype &x , Statetype &dxdt , const double )

The first argument should be const.

I'm not sure this is the problem though - the error message is rather cryptic...

[jannisehrlich]

Making the first argument of the operator() const doesnot help. Indeed the error message is rather cryptic - that's why I don't get it resolved...

[ddemidov]

It looks like the real problem is that there is no algebra dispatcher for boost::multi_array<T,N>. So you either have to specify the algebra explicitly, as in

#include <iostream>
#include <complex>
#include <cmath>
#include <boost/multi_array.hpp>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/algebra/multi_array_algebra.hpp>

typedef boost::multi_array< double , 2 > Statetype;

class harm_osc {
private:
    double m_gam;

public:
    harm_osc( double gam ) : m_gam(gam) { }

    void operator() ( const Statetype &x , Statetype &dxdt , const double )
    {
        dxdt[0][0] = x[1][0];
        dxdt[1][0] = -x[0][0] - m_gam*x[1][0];
    }
};

int main(){
    Statetype x;
    x.resize(boost::extents[2][1]);
    x[0][0]=1.0;
    x[1][0]=0.0;

    boost::numeric::odeint::runge_kutta4 <
        Statetype, double,
        Statetype, double,
        boost::numeric::odeint::multi_array_algebra
        > stepper;

    harm_osc ho(0.15);
    double dt = 0.1;

    for( double t=0.0 ; t<10.0 ; t+= dt ){
        stepper.do_step( ho , x , t , dt );
        std::cout << t << '\t' << x[0][0] << '\t' << x[1][0] << '\n';
    }
}

or provide the dispatcher yourself:

#include <iostream>
#include <complex>
#include <cmath>
#include <boost/multi_array.hpp>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/algebra/multi_array_algebra.hpp>

typedef boost::multi_array< double , 2 > Statetype;

namespace boost { namespace numeric { namespace odeint {
template< class T , size_t N >
struct algebra_dispatcher< boost::multi_array< T , N > >
{
    typedef multi_array_algebra algebra_type;
};
} } }

class harm_osc {
private:
    double m_gam;

public:
    harm_osc( double gam ) : m_gam(gam) { }

    void operator() ( const Statetype &x , Statetype &dxdt , const double )
    {
        dxdt[0][0] = x[1][0];
        dxdt[1][0] = -x[0][0] - m_gam*x[1][0];
    }
};

int main(){
    Statetype x;
    x.resize(boost::extents[2][1]);
    x[0][0]=1.0;
    x[1][0]=0.0;

    boost::numeric::odeint::runge_kutta4 <Statetype> stepper;

    harm_osc ho(0.15);
    double dt = 0.1;

    for( double t=0.0 ; t<10.0 ; t+= dt ){
        stepper.do_step( ho , x , t , dt );
        std::cout << t << '\t' << x[0][0] << '\t' << x[1][0] << '\n';
    }
}

[ddemidov]

206 should fix this.

[jannisehrlich]

Thanks that is working. But for the case that I have the multi_array inside of a class - as in the other part of my code - I have to write my own algebra dispatcher?

[ddemidov]

My guess is you either have to provide your own algebra, or make your state type compatible with one of the provided algebras.

[jannisehrlich]

Ok then I will change and use the provided algebras. Thanks for your help