Blaze Vectors as state type for odeint library with adaptive integration

[PhilSche]

Hi Klaus,

I am trying to use Blaze vectors as a state type for the integration with odeint. The integration with constant step size works perfectly fine, but once I switch to the adaptive routine my code does not compile anymore due to type issues which I am unable to track down.

You can find a minimal example which demonstrates the issue below:

#include <blaze/Math.h>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/blaze/blaze_algebra_dispatcher.hpp>
#include <boost/numeric/odeint/external/blaze/blaze_resize.hpp>

using namespace blaze;
using namespace boost::numeric::odeint;


typedef blaze::DynamicVector<double> State_Type;


class ODEsystem {
    public:
        ODEsystem() { };    
        void operator() (const State_Type &y, State_Type &dydt, const double /*t*/){ 
            const State_Type B{0., 0., 1.};
            dydt = cross(y, B);
        };
};



int main(int argc, char** argv) {
    State_Type y{1., 0., 0.};   // initial condition
    ODEsystem sys;

    typedef runge_kutta_dopri5< State_Type, double, State_Type, double, vector_space_algebra > stepper_type;
    typedef controlled_runge_kutta< stepper_type > controlled_stepper_type;
    controlled_runge_kutta< stepper_type > controlled_stepper = make_controlled< stepper_type >(1e-6, 1e-6);    // stepper with adaptive step size control
    runge_kutta_dopri5< State_Type, double, State_Type, double, vector_space_algebra > stepper;     // stepper used for constant step size

    integrate_const(stepper, sys, y, 0., 100., 1e-3);   // working

    integrate_adaptive(controlled_stepper, sys, y, 0., 100., 1e-3); // not working

    return 0;
}

Does the issue appear due to a missing compatibility between both libraries or is it simple me using the odeint library in a wrong way?

I'm compiling with g++ 8.2 using Blaze 3.4 and Boost 1.68.

Best regards, Philipp

[Klaus Iglberger]

Hi Philipp!

The above code compiles with the following error message:

/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:89:24: error: no matching member function for call to 'norm_inf'
        return algebra.norm_inf( x_err );

This message indeed indicates an incompatibility between the two libraries or a misuse of odeint. Perhaps following the individual steps in the odeint harmonic oscillator example help you to identify the problem. From our point of view, however, this is definitely not a bug in the Blaze library and we will close the issue.

Best regards,

Klaus!

[PhilSche]

Hi Klaus,

I managed to fix the the issue by including the following lines in the code:

state_type operator+( const state_type &lhs , const double rhs ) {
    state_type new_state = lhs;
    for (auto i = 0; i < lhs.size(); ++i){
        new_state[i] += rhs;
    }
    return new_state;
}

state_type operator+( const double lhs , const state_type &rhs ) {
    state_type new_state = rhs;
    for (auto i = 0; i < rhs.size(); ++i){
        new_state[i] += lhs;
    }
    return new_state;
}


namespace boost {
    namespace numeric {
        namespace odeint {
            template<>
            struct vector_space_norm_inf<state_type> {
                typedef double result_type;
                auto operator()(const state_type &u) const {
                    return blaze::max(u);
                }
            };
}
}
}

It appears that blaze does not define a behavior for vector/scalar addition, which seems to be required by odeint.

I assume that my implementation is very inefficient. Is it possible to include this functionality it in one of the next releases?

Best regards, Philipp

[Klaus Iglberger]

Hi Philipp!

You are correct, Blaze currently does provide a vector/scalar addition. In order to provide us a reminder to add the feature it would be great if you could create a "proposal" issue. Please provide a short description of the feature and possibly some example code (the one you posted above would be sufficient). Thanks,

Best regards,

Klaus!