Skip to content

Instantly share code, notes, and snippets.

@ddemidov
Last active May 29, 2017 19:22
Show Gist options
  • Save ddemidov/8331121 to your computer and use it in GitHub Desktop.
Save ddemidov/8331121 to your computer and use it in GitHub Desktop.
odeint with boost.compute backend
#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;
}
}
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment