Параметры шаблона `boost::numeric::odeint::runge_kutta-X` совместимы с CUDA/OpenMP

Question

Подпись типа для класса stepper, который я использую, суммирована здесь:

http://www.boost.org/doc/libs/1_56_0/libs/numeric/odeint/doc/html/boost/numeric/odeint/runge_kutta_dopri5.html

Это может быть реализовано следующим образом:

 boost::numeric::odeint::runge_kutta_dopri5< state_type_ > stepper;

Все идет нормально. Оно работает.

Я планирую перенести мою программу в Cuda (используя Thrust), а затем в OpenMP. Я изменил объявление на следующее:

boost::numeric::odeint::runge_kutta_dopri5< state_type_
        , double
        , state_type_
        , double
        , boost::numeric::odeint::vector_space_algebra 
        > stepper;

Я следовал за решением этой проблемы, но это не компилируется.

In file included from /usr/include/boost/numeric/odeint/stepper/euler.hpp:26:
/usr/include/boost/numeric/odeint/algebra/default_operations.hpp:87:27: error: invalid operands to binary expression ('double' and 'const std::vector<double, std::allocator<double> >')
            t1 = m_alpha1 * t2 + m_alpha2 * t3;
                 ~~~~~~~~ ^ ~~

Мне интересно, какой самый портативный способ объявить степпер, чтобы минимальные изменения потребовались позже при портировании на cuda.

Answer 1

Это зависит от того, что вы хотите сделать. Если вы хотите использовать Thrust, вам нужно изменить объявление на

boost::numeric::odeint::runge_kutta_dopri5<
    state_type , double , state_type , double ,
    thrust_algebra , thrust_operations >;

thrust_algebra и thrust_operations убедитесь, что все вычисления перенаправлены на соответствующие thrust::for_each вызовы, где используются заархивированные итераторы. Если вы хотите использовать высокоуровневую библиотеку линейной алгебры высокого уровня, которая работает на графическом процессоре (например, VexCL или ViennaCL), вы можете использовать вышеуказанное объявление и изменить только state_type к правильному типу, например vexcl::vector< double >, vector_space_algebra предполагает, что ваш state_type может обрабатывать такие операции, как y = a1*x1 + a2*x2, что имеет место для VexCL и ViennaCL из-за использования шаблонов выражений. Вы также можете посмотреть здесь.