odeint with boost.compute backend

lorenz.cpp

#include <iostream>
#include <iomanip>

#include <boost/timer/timer.hpp>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/compute/compute.hpp>
#include <boost/compute.hpp>

namespace odeint = boost::numeric::odeint;
namespace bc     = boost::compute;

typedef bc::vector<double> state_type;

//---------------------------------------------------------------------------
struct lorenz_system {
    size_t n;
    double sigma, b;
    const state_type &R;

    lorenz_system(size_t n, const state_type &R,
            double sigma = 10.0, double b = 8.0 / 3.0
            ) : n(n), R(R), sigma(sigma), b(b)
    {}

    template <class State, class Deriv>
    void operator()(const State &x, Deriv &dxdt, double t) const {
        using bc::lambda::make_tuple;

        auto X  = bc::lambda::get<0>(bc::_1);
        auto Y  = bc::lambda::get<1>(bc::_1);
        auto Z  = bc::lambda::get<2>(bc::_1);
        auto r  = bc::lambda::get<3>(bc::_1);
        auto dX = bc::lambda::get<4>(bc::_1);
        auto dY = bc::lambda::get<5>(bc::_1);
        auto dZ = bc::lambda::get<6>(bc::_1);

        bc::for_each(
                bc::make_zip_iterator(
                    boost::make_tuple(
                        x.begin(),
                        x.begin() + n,
                        x.begin() + 2 * n,
                        R.begin(),
                        dxdt.begin(),
                        dxdt.begin() + n,
                        dxdt.begin() + 2 * n
                        )
                    ),
                bc::make_zip_iterator(
                    boost::make_tuple(
                        x.begin() + n,
                        x.begin() + 2 * n,
                        x.end(),
                        R.end(),
                        dxdt.begin() + n,
                        dxdt.begin() + 2 * n,
                        dxdt.end()
                        )
                    ),
                make_tuple(
                        dX = sigma * (Y - X),
                        dY = r * X - Y - X * Z,
                        dZ = -b * Z + X * Y
                        )
                );
    }
};

//---------------------------------------------------------------------------
int main(int argc, char *argv[]) {

    const size_t n = argc > 1 ? atoi(argv[1]) : 1024;
    const double dt = 0.01;
    const double t_max = 10.0;

    std::cout << bc::system::default_device().name() << std::endl;

    std::vector<double> r(n);
    const double Rmin = 0.1, Rmax = 50.0;
    for(size_t i = 0; i < n; ++i)
        r[i] = Rmin + i * (Rmax - Rmin) / (n - 1);

    state_type R(n);
    state_type X(3 * n);

    bc::copy(r.begin(), r.end(), R.begin());

    // initialize x,y,z
    bc::fill(X.begin(), X.end(), 10.0);

    typedef odeint::runge_kutta4_classic< state_type > stepper_type;

    lorenz_system lorenz(n, R);

    boost::timer::cpu_timer timer;

    odeint::integrate_const(stepper_type(), boost::ref(lorenz), X, 0.0, t_max, dt);

    bc::system::finish();
    boost::timer::cpu_times t = timer.elapsed();
    std::cout << boost::timer::format(t, 8, "%w") << std::endl;

    if (n <= 1024) {
        std::vector<double> x(3 * n);
        bc::copy(X.begin(), X.end(), x.begin());

        std::cout << "x = {" << std::setprecision(6);
        for(size_t i = 0; i < n; ++i) {
            if (i % 2 == 0) std::cout << "\n" << std::setw(6) << i << ":";
            std::cout << std::scientific << " (";
            for(size_t j = 0; j < 3; ++j)
                std::cout << std::setw(14) << x[j * n + i];
            std::cout << ")";
        }
        std::cout << "\n}" << std::endl;
    }
}

Makefile

lorenz: lorenz.cpp
	g++ -o $@ $^ -std=c++11 -O3 \
	    -I${HOME}/work/odeint-v2 \
	    -I${HOME}/work/opencl/compute/include/ \
	    -lOpenCL -lboost_system -lboost_timer