В предыдущем разделе показано использование odeint в сочетании с Thrust. В этом разделе мы покажем, как можно использовать OpenCL с odeint. Суть odeint заключается не в реализации собственных низкоуровневых структур данных и алгоритмов, а в использовании высокоуровневых библиотек, выполняющих эту задачу. Здесь мы будем использовать фреймворк VexCL для использования OpenCL. VexCL - хорошая библиотека для общих вычислений, и она использует шаблоны с высокой степенью экспрессии. С помощью VexCL можно писать очень компактное и выразительное приложение.

	Note
	vexcl нуждается в C++11! Таким образом, вы должны скомпилировать с поддержкой C++11.

Для использования VexCL необходимо включить один дополнительный заголовок, который включает в себя типы данных и алгоритмы от vexcl и адаптацию к odeint. Адаптация к odeint означает здесь только адаптацию функционала изменения размера VexCL к odeint.

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

Для демонстрации использования VexCL мы интегрируем ансамбль системы Лоренца. Пример очень похож на изучение параметров системы Лоренца в предыдущем разделе, за исключением того, что мы не вычисляем показатели Ляпунова. Опять же, мы варьируем параметр R системы Лоренца, параллельно решая целый ансамбль систем Лоренца (каждая с другим параметром R). Во-первых, мы определяем тип состояния и тип вектора

typedef vex::vector< double >    vector_type;
typedef vex::multivector< double, 3 > state_type;

vector_type используется для представления параметра R. state_type является мультивектором из трех подвекторов и используется для представления. Первый компонент этого мультивектора представляет все компоненты системы Лоренца x, а второй - все компоненты системы Лоренца y и третий - все компоненты z. Компоненты этого вектора могут быть получены через

auto &x = X(0);
auto &y = X(1);
auto &z = X(2);

Как уже упоминалось, VexCL поддерживает шаблоны выражения, и мы будем использовать их для реализации функции системы для ансамбля Лоренца:

const double sigma = 10.0;
const double b = 8.0 / 3.0;
struct sys_func
{
    const vector_type &R;
    sys_func( const vector_type &_R ) : R( _R ) { }
    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {
        dxdt(0) = -sigma * ( x(0) - x(1) );
        dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
        dxdt(2) = - b * x(2) + x(0) * x(1);
    }
};

Это очень просто, не так ли? Эти три маленькие строки делают все вычисления для вас. Нет необходимости писать собственные ядра OpenCL. VexCL делает все для вас. Далее нужно написать основную заявку. Мы инициализируем вектор параметров (R) и начальное состояние. Обратите внимание, что VexCL требует vector_space_algebra, но это автоматически выводится и настраивается одеином внутри, поэтому при инстанцировании степпера мы должны указать только state_type:

// setup the opencl context
vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
std::cout << ctx << std::endl;
// set up number of system, time step and integration time
const size_t n = 1024 * 1024;
const double dt = 0.01;
const double t_max = 1000.0;
// initialize R
double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
std::vector<double> x( n * 3 ) , r( n );
for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
vector_type R( ctx.queue() , r );
// initialize the state of the lorenz ensemble
state_type X(ctx.queue(), n);
X(0) = 10.0;
X(1) = 10.0;
X(2) = 10.0;
// create a stepper
runge_kutta4< state_type > stepper;
// solve the system
integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );