Открытый одеинт Поддержка MP реализована как внешний бэкэнд, который должен быть включен вручную. В зависимости от компилятора могут потребоваться дополнительные флаги, т.е.<-fopenmp>для GCC.

#include <omp.h>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/openmp/openmp.hpp>

В самом простом подходе к параллелизации с OpenMP мы используем стандарт<vector>в качестве типа состояния:

typedef std::vector< double > state_type;

Мы инициализируем состояние с некоторыми случайными данными:

size_t N = 131101;
state_type x( N );
boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
boost::random::mt19937 engine( 0 );
generate( x.begin() , x.end() , boost::bind( distribution , engine ) );

Теперь нам нужно настроить степпер для использования бэкэнда OpenMP. Это делается путем явного предоставления<openmp_range_algebra>в качестве параметра шаблона для степпера. Эта алгебра требует, чтобы тип состояния был моделью диапазона случайного доступа и использовался алгеброй из нескольких потоков.

typedef runge_kutta4<
                  state_type , double ,
                  state_type , double ,
                  openmp_range_algebra
                > stepper_type;

В дополнение к предоставлению степпера с OpenMP-параллелизацией нам также нужна функция параллельной системы для использования доступных ядер. Здесь показана простая одномерная цепь фазовых осцилляторов с соединением ближайшего соседа:

struct phase_chain
{
    phase_chain( double gamma = 0.5 )
    : m_gamma( gamma ) { }
    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
    {
        const size_t N = x.size();
        #pragma omp parallel for schedule(runtime)
        for(size_t i = 1 ; i < N - 1 ; ++i)
        {
            dxdt[i] = coupling_func( x[i+1] - x[i] ) +
                      coupling_func( x[i-1] - x[i] );
        }
        dxdt[0  ] = coupling_func( x[1  ] - x[0  ] );
        dxdt[N-1] = coupling_func( x[N-2] - x[N-1] );
    }
    double coupling_func( double x ) const
    {
        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
    }
    double m_gamma;
};

	Note
	В бэкэндах OpenMP функция системы всегда будет называться последовательно из потока, используемого для запуска интеграции.

Наконец, мы выполняем интеграцию, используя одну из функций интеграции от odeint. Как видите, параллелизация полностью скрыта в степпере и функции системы. Открыть Депутат позаботится о распределении работы между нитями и присоединится к ним автоматически.

integrate_n_steps( stepper_type() , phase_chain( 1.2 ) ,
                   x , 0.0 , 0.01 , 100 );

После интеграции данные могут быть немедленно доступны и обработаны. Обратите внимание, что вы можете указать расписание OpenMP, позвонив<omp_set_schedule>в начале программы:

int chunk_size = N/omp_get_max_threads();
omp_set_schedule( omp_sched_static , chunk_size );

См.openmp/phase_chain.cppДля полного примера.

Split state

Для продвинутых случаев одеинт предлагает другой подход к использованию OpenMP, который позволяет более точно контролировать параллелизацию. Например, для данных нечетного размера, где границы потоков OpenMP не совпадают с линиями кэша и повреждают производительность, может быть целесообразно скопировать данные из непрерывного<vector<T>>в отдельные, индивидуально выровненные векторы. Для этого odeint предоставляет тип<openmp_state<T>>, по существу псевдоним для<vector<vector<T>>>.

Здесь инициализация выполняется с<vector<double>>, но затем мы используем функцию odeint<split>для заполнения<openmp_state>. Расщепление выполняется таким образом, что размеры отдельных областей отличаются максимум на 1, чтобы сделать вычисления как можно более однородными.

const size_t N = 131101;
vector<double> x( N );
boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
boost::random::mt19937 engine( 0 );
generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
const size_t blocks = omp_get_max_threads();
state_type x_split( blocks );
split( x , x_split );

Конечно, нужно изменить функцию системы, чтобы справиться с<openmp_state>. Отметим, что каждый субрегион государства вычисляется в едином задании, но на границах считывается доступ к соседним регионам.

struct phase_chain_omp_state
{
    phase_chain_omp_state( double gamma = 0.5 )
    : m_gamma( gamma ) { }
    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
    {
        const size_t N = x.size();
        #pragma omp parallel for schedule(runtime)
        for(size_t n = 0 ; n < N ; ++n)
        {
            const size_t M = x[n].size();
            for(size_t m = 1 ; m < M-1 ; ++m)
            {
                dxdt[n][m] = coupling_func( x[n][m+1] - x[n][m] ) +
                             coupling_func( x[n][m-1] - x[n][m] );
            }
            dxdt[n][0] = coupling_func( x[n][1] - x[n][0] );
            if( n > 0 )
            {
                dxdt[n][0] += coupling_func( x[n-1].back() - x[n].front() );
            }
            dxdt[n][M-1] = coupling_func( x[n][M-2] - x[n][M-1] );
            if( n < N-1 )
            {
                dxdt[n][M-1] += coupling_func( x[n+1].front() - x[n].back() );
            }
        }
    }
    double coupling_func( double x ) const
    {
        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
    }
    double m_gamma;
};

Используя<openmp_state<T>>тип состояния автоматически выбирает<openmp_algebra>, который выполняет внутренние вычисления odeint на параллельных областях. Следовательно, ручная конфигурация степпера не требуется. В конце интеграции мы используем<unsplit>для объединения субрегионов в единый вектор.

integrate_n_steps( runge_kutta4<state_type>() , phase_chain_omp_state( 1.2 ) ,
                   x_split , 0.0 , 0.01 , 100 );
unsplit( x_split , x );

	Note
	На самом деле вам не нужно использовать<openmp_state<T>>для расширенных вариантов использования,<openmp_algebra>является просто псевдонимом для<openmp_nested_algebra<range_algebra>>и поддерживает любую модель диапазона случайного доступа в качестве внешнего типа параллельного состояния и будет использовать данную алгебру на своих элементах.

См.openmp/phase_chain_omp_state.cppдля полного примера.

Для расширения параллельных вычислений на нескольких машинах мы можем использовать MPI.

Реализация системной функции аналогична варианту OpenMP с разделенными данными, основное отличие заключается в том, что, хотя OpenMP использует модель нереста / соединения, где все явно не параллельные выполняется только в основном потоке, в модели MPI каждый узел входит в метод<main()>независимо, расходясь на основе его ранга и синхронизируя через пропускание сообщений и явные барьеры.

Поддержка MPI odeint также реализована как внешний бэкэнд. В зависимости от реализации MPI код может быть скомпилирован с помощью<mpic++>.

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

Вместо того, чтобы читать данные другого потока, мы асинхронно отправляем и получаем соответствующие данные из соседних узлов, выполняя некоторые вычисления в промежуточный период, чтобы скрыть задержку.

struct phase_chain_mpi_state
{
    phase_chain_mpi_state( double gamma = 0.5 )
    : m_gamma( gamma ) { }
    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
    {
        const size_t M = x().size();
        const bool have_left = x.world.rank() > 0,
                   have_right = x.world.rank() < x.world.size()-1;
        double x_left, x_right;
        boost::mpi::request r_left, r_right;
        if( have_left )
        {
            x.world.isend( x.world.rank()-1 , 0 , x().front() ); // send to x_right
            r_left = x.world.irecv( x.world.rank()-1 , 0 , x_left ); // receive from x().back()
        }
        if( have_right )
        {
            x.world.isend( x.world.rank()+1 , 0 , x().back() ); // send to x_left
            r_right = x.world.irecv( x.world.rank()+1 , 0 , x_right ); // receive from x().front()
        }
        for(size_t m = 1 ; m < M-1 ; ++m)
        {
            dxdt()[m] = coupling_func( x()[m+1] - x()[m] ) +
                        coupling_func( x()[m-1] - x()[m] );
        }
        dxdt()[0] = coupling_func( x()[1] - x()[0] );
        if( have_left )
        {
            r_left.wait();
            dxdt()[0] += coupling_func( x_left - x().front() );
        }
        dxdt()[M-1] = coupling_func( x()[M-2] - x()[M-1] );
        if( have_right )
        {
            r_right.wait();
            dxdt()[M-1] += coupling_func( x_right - x().back() );
        }
    }
    double coupling_func( double x ) const
    {
        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
    }
    double m_gamma;
};

Аналогично<openmp_state<T>>мы используем<mpi_state<InnerState<T>>>, который автоматически выбирает<mpi_nested_algebra>и соответствующую MPI-незаметную внутреннюю алгебру (поскольку наше внутреннее состояние является<vector>, внутренняя алгебра будет<range_algebra>, как в примере OpenMP).

typedef mpi_state< vector<double> > state_type;

В основной программе мы строим<communicator>, который сообщает нам<size>кластера и текущий узел<rank>в нем. Мы генерируем входные данные только на главном узле, избегая ненужной работы на других узлах. Вместо того, чтобы просто копировать куски,<split>действует здесь как коллективная функция MPI и отправляет / получает области от хозяина к каждому рабу. Входной аргумент игнорируется на рабах, но главный узел получает область в своем выходе и будет участвовать в вычислении.

boost::mpi::environment env( argc , argv );
boost::mpi::communicator world;
const size_t N = 131101;
vector<double> x;
if( world.rank() == 0 )
{
    x.resize( N );
    boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
    boost::random::mt19937 engine( 0 );
    generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
}
state_type x_split( world );
split( x , x_split );

Теперь, когда<x_split>содержит (только) локальный фрагмент для каждого узла, мы начинаем интеграцию.

Чтобы распечатать результат на главном узле, мы отправляем обработанные данные обратно с помощью<unsplit>.

integrate_n_steps( runge_kutta4<state_type>() , phase_chain_mpi_state( 1.2 ) ,
                   x_split , 0.0 , 0.01 , 100 );
unsplit( x_split , x );

Note

<mpi_nested_algebra::for_each><N>не использует никаких MPI-конструкций, он просто вызывает внутреннюю алгебру на локальном участке, и функция системы также не защищена никакими барьерами, поэтому, если вы вручную не разместите какие-либо (например, в случаях исследования параметров, когда элементы полностью независимы), вы можете увидеть, что узлы расходятся, возвращаясь от этого вызова в разное время.

См.mpi/phase_chain.cppдля полного примера.