|
|
|
Using CUDA (or OpenMP, TBB, ...) via ThrustBoost , Chapter 1. Boost.Numeric.Odeint , Tutorial
|
|||||||||||||||||||||
 | Home | Libraries | People | FAQ | More |
Современные графические карты (графические процессоры - графические процессоры) могут использоваться для ускорения производительности алгоритмов, затрачивающих время, посредством массивной параллелизации. Они предназначены для выполнения многих операций параллельно. odeint может использовать мощность графических процессоров с помощью CUDA иThrust, который является STL-подобным интерфейсом для собственного API CUDA.
![[Important]](/img/important.png) | Important |
|---|---|
Thrust также поддерживает параллелизацию с использованием блоков OpenMP и Intel Threading Building Blocks (TBB). Вы можете переключаться между CUDA, OpenMP и TBB параллелизациями с помощью простого коммутатора компилятора. Следовательно, это также обеспечивает простой способ получить базовую параллелизацию OpenMP в odeint. Приведенные ниже примеры сосредоточены на параллелизации GPU. |
Для использования odeint с CUDA необходимо учитывать несколько моментов. Прежде всего, проблему нужно правильно выбрать. Совершенно бессмысленно пытаться распараллеливать код для трехмерной системы, он просто слишком мал и не стоит усилий. Один вызов функции (выполнение ядра) на GPU происходит медленно, но вы можете выполнить операцию на огромном наборе данных только одним вызовом. Мы испытали, что размер вектора, над которым проводится параллелизация, должен быть порядка106, чтобы в полной мере использовать GPU. Во-вторых, вы должны использоватьАлгоритмы и функторы Thrustпри реализации Rhs ODE. Это может быть сложно, поскольку включает в себя некоторые знания функционального программирования.
Типичными приложениями для CUDA и odeint являются большие системы, такие как решетки или дискретизации ФДЭ и исследования параметров. В настоящее время мы представляем три примера, которые показывают, как мощность графических процессоров может использоваться в сочетании с odeint.
| Important |
|---|---|
Полная мощность CUDA доступна только для действительно больших систем, где число связанных обычных дифференциальных уравнений имеет порядокN=106или больше. Для небольших систем процессор обычно намного быстрее. Вы также можете интегрировать ансамбль различных несвязанных ODE параллельно, как показано в последнем примере. |
Первым примером является ансамбль фазовых осцилляторов из предыдущего раздела:
dφk / dt = ωk + ε / N Σj sin( φj - φk ).
Он имеет фазовый переход наи #949; = 2в пределе бесконечного числа осцилляторовN. В случае конечногоNэтот переход размазан, но все еще хорошо виден.
Тягаи CUDA идеально подходят для таких задач, когда требуется большое количество частиц (осцилляторов). Начнем с определения типа государства<thrust::device_vector>. Содержание этого вектора живет на GPU. Если вы не знакомы с этим, мы рекомендуем прочитать.НачалоразделаСсылка на сайт.
//change this to float if your device does not support double computation typedef double value_type; //change this to host_vector< ... > of you want to run on CPU typedef thrust::device_vector< value_type > state_type; // typedef thrust::host_vector< value_type > state_type;
Thrust следует подходу функционального программирования. Если вы хотите выполнить вычисление на GPU, вам обычно приходится вызывать глобальную функцию, такую как<thrust::for_each>,<thrust::reduce>, ... с соответствующим локальным функтором, который выполняет основную операцию. Примером является
struct add_two { template< class T > __host__ __device__ void operator()( T &t ) const { t += T( 2 ); } }; // ... thrust::for_each( x.begin() , x.end() , add_two() );
Этот код обычно добавляет два к каждому элементу в контейнере<x>.
Для интеграции фазового осциллятора необходим ансамбль
The mean field is calculated in a class mean_field_calculator
struct mean_field_calculator { struct sin_functor : public thrust::unary_function< value_type , value_type > { __host__ __device__ value_type operator()( value_type x) const { return sin( x ); } }; struct cos_functor : public thrust::unary_function< value_type , value_type > { __host__ __device__ value_type operator()( value_type x) const { return cos( x ); } }; static std::pair< value_type , value_type > get_mean( const state_type &x ) { value_type sin_sum = thrust::reduce( thrust::make_transform_iterator( x.begin() , sin_functor() ) , thrust::make_transform_iterator( x.end() , sin_functor() ) ); value_type cos_sum = thrust::reduce( thrust::make_transform_iterator( x.begin() , cos_functor() ) , thrust::make_transform_iterator( x.end() , cos_functor() ) ); cos_sum /= value_type( x.size() ); sin_sum /= value_type( x.size() ); value_type K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum ); value_type Theta = atan2( sin_sum , cos_sum ); return std::make_pair( K , Theta ); } };
Внутри этого класса определены две членские структуры<sin_functor>и<cos_functor>. Они вычисляют синус и косинус значения и используются в итераторе преобразования для вычисления суммыгреха [φ& #8203;k]иcos [φ& #8203;k]. Классификаторы<__host__>и<__device__>являются CUDA-специфичными и определяют функцию или оператора, которые могут выполняться как на GPU, так и на CPU. Линия
value_type sin_sum = thrust::reduce( thrust::make_transform_iterator( x.begin() , sin_functor() ) , thrust::make_transform_iterator( x.end() , sin_functor() ) );
выполняет вычисление этой синусовой суммы на GPU (или на CPU, в зависимости от конфигурации тяги).
Функция системы определяется через
class phase_oscillator_ensemble { public: struct sys_functor { value_type m_K , m_Theta , m_epsilon; sys_functor( value_type K , value_type Theta , value_type epsilon ) : m_K( K ) , m_Theta( Theta ) , m_epsilon( epsilon ) { } template< class Tuple > __host__ __device__ void operator()( Tuple t ) { thrust::get<2>(t) = thrust::get<1>(t) + m_epsilon * m_K * sin( m_Theta - thrust::get<0>(t) ); } }; // ... void operator() ( const state_type &x , state_type &dxdt , const value_type dt ) const { std::pair< value_type , value_type > mean_field = mean_field_calculator::get_mean( x ); thrust::for_each( thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_omega.begin() , dxdt.begin() ) ), thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_omega.end() , dxdt.end()) ) , sys_functor( mean_field.first , mean_field.second , m_epsilon ) ); } // ... };
Этот класс используется в методах<do_step>и<integrate>. Он определяет структуру члена<sys_functor>для r.h.s. каждого отдельного осциллятора и<operator()>для использования в степперах и интеграторах odeint. Функтор вычисляет сначала среднее поле& #966;& #8203;kи во-вторых вычисляет все r.h.s. ODE с использованием этого среднего поля. Как хорошо<thrust::tuple>и<thrust::zip_iterator>играть вместе.
Теперь мы готовы собрать все воедино. Все, что нам нужно сделать, чтобы подготовить одеинт к использованию графического процессора, это параметризировать степпер с<state_type>и<value_type>:
typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;
![[Note]](/img/note.png) | Note |
|---|---|
Мы специально определили четыре параметра шаблона, потому что мы должны переопределить значение параметра по умолчанию< |
Вы также можете использовать контролируемый или плотный выходной степпер, например.
typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
Затем легко интегрировать фазовый ансамбль, создав экземпляр класса rhs и используя функцию интеграции:
phase_oscillator_ensemble ensemble( N , 1.0 );
size_t steps1 = integrate_const( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
Мы должны использовать<boost::ref>здесь, чтобы передать класс rhs в качестве эталона, а не по стоимости. Это гарантирует, что естественные частоты каждого осциллятора не копируются при вызове<integrate_const>. В полном примере сравниваются производительность и результаты решения Runge-Kutta-4 и Dopri5.
Полный пример можно найти наphase_oscillator_example.cu.
Следующим примером является большая одномерная цепочка фазовых осцилляторов, связанных с ближайшим соседом, со следующими уравнениями движения:
d φk / dt = ωk + sin( φk+1 - φk ) + sin( φk - φk-1)
В принципе, мы можем использовать все методы из предыдущего примера ансамбля фазовых осцилляторов, но мы должны уделять особое внимание соединению осцилляторов. Чтобы эффективно реализовать соединение, вы можете использовать очень элегантный способ, используя итератор перестановки Thrust. Пермутационный итератор ведет себя как обычный итератор на векторе, но он не повторяется в соответствии с обычным порядком элементов. Он скорее повторяется вдоль некоторой перестановки элементов, определенных некоторой картой индекса. Чтобы реализовать связь с ближайшим соседом, мы создаем один итератор перестановок, который проходит один шаг позади обычного итератора, и другой итератор перестановок, который проходит один шаг впереди. Полный класс системы:
//change this to host_vector< ... > if you want to run on CPU typedef thrust::device_vector< value_type > state_type; typedef thrust::device_vector< size_t > index_vector_type; //typedef thrust::host_vector< value_type > state_type; //typedef thrust::host_vector< size_t > index_vector_type; class phase_oscillators { public: struct sys_functor { template< class Tuple > __host__ __device__ void operator()( Tuple t ) // this functor works on tuples of values { // first, unpack the tuple into value, neighbors and omega const value_type phi = thrust::get<0>(t); const value_type phi_left = thrust::get<1>(t); // left neighbor const value_type phi_right = thrust::get<2>(t); // right neighbor const value_type omega = thrust::get<3>(t); // the dynamical equation thrust::get<4>(t) = omega + sin( phi_right - phi ) + sin( phi - phi_left ); } }; phase_oscillators( const state_type &omega ) : m_omega( omega ) , m_N( omega.size() ) , m_prev( omega.size() ) , m_next( omega.size() ) { // build indices pointing to left and right neighbours thrust::counting_iterator<size_t> c( 0 ); thrust::copy( c , c+m_N-1 , m_prev.begin()+1 ); m_prev[0] = 0; // m_prev = { 0 , 0 , 1 , 2 , 3 , ... , N-1 } thrust::copy( c+1 , c+m_N , m_next.begin() ); m_next[m_N-1] = m_N-1; // m_next = { 1 , 2 , 3 , ... , N-1 , N-1 } } void operator() ( const state_type &x , state_type &dxdt , const value_type dt ) { thrust::for_each( thrust::make_zip_iterator( thrust::make_tuple( x.begin() , thrust::make_permutation_iterator( x.begin() , m_prev.begin() ) , thrust::make_permutation_iterator( x.begin() , m_next.begin() ) , m_omega.begin() , dxdt.begin() ) ), thrust::make_zip_iterator( thrust::make_tuple( x.end() , thrust::make_permutation_iterator( x.begin() , m_prev.end() ) , thrust::make_permutation_iterator( x.begin() , m_next.end() ) , m_omega.end() , dxdt.end()) ) , sys_functor() ); } private: const state_type &m_omega; const size_t m_N; index_vector_type m_prev; index_vector_type m_next; };
Обратите внимание, насколько легко можно получить значение для левого и правого соседнего осциллятора в системном функторе с помощью итераторов перестановок. Но вызов функции<thrust::for_each>выглядит относительно сложным. Каждый термин r.h.s. ODE напоминает один итератор, упакованный точно так же, как он упакован в функторе выше.
Теперь мы собрали все вместе. Мы создаем случайные начальные условия и уменьшаем частоты так, чтобы получить синхронизацию. Мы копируем частоты и начальные условия на устройство и, наконец, инициализируем и выполняем интеграцию. В результате мы просто записываем текущее состояние, следовательно, фазу каждого осциллятора.
// create initial conditions and omegas on host: vector< value_type > x_host( N ); vector< value_type > omega_host( N ); for( size_t i=0 ; i<N ; ++i ) { x_host[i] = 2.0 * pi * drand48(); omega_host[i] = ( N - i ) * epsilon; // decreasing frequencies } // copy to device state_type x = x_host; state_type omega = omega_host; // create stepper runge_kutta4< state_type , value_type , state_type , value_type > stepper; // create phase oscillator system function phase_oscillators sys( omega ); // integrate integrate_const( stepper , sys , x , 0.0 , 10.0 , dt ); thrust::copy( x.begin() , x.end() , std::ostream_iterator< value_type >( std::cout , "\n" ) ); std::cout << std::endl;
Полный пример можно найти наphase_oscillator_chain.cu.
Другим важным примером использованияThrustи CUDA являются исследования параметров относительно небольших систем. Рассмотрим, например, трехмерную систему Лоренца из примера хаотических систем в предыдущем разделе, которая имеет три параметра. Если вы хотите изучить поведение этой системы по различным параметрам, вам обычно приходится интегрировать систему для многих значений параметров. Используя тягу и одеинт, вы можете сделать эту интеграцию параллельно, поэтому вы интегрируете целый ансамбль систем Лоренца, где каждая индивидуальная реализация имеет разное значение параметра.
Ниже мы покажем, как можно использоватьThrustдля интеграции вышеупомянутого ансамбля систем Лоренца. Мы будем изменять только параметри #946;, но легко изменить другие параметры или даже два или все три параметра. Кроме того, мы будем использовать самый большой показатель Ляпунова для количественной оценки поведения системы.
Мы начинаем с определения диапазона параметров, которые мы хотим изучить. Состояние_тип снова является<thrust::device_vector<value_type
>>.
vector< value_type > beta_host( N ); const value_type beta_min = 0.0 , beta_max = 56.0; for( size_t i=0 ; i<N ; ++i ) beta_host[i] = beta_min + value_type( i ) * ( beta_max - beta_min ) / value_type( N - 1 ); state_type beta = beta_host;
Следующее, что мы должны реализовать, — это система Лоренца без возмущений. Позже также реализуется система с возмущениями для вычисления Ляпуновского экспонента. Мы будем использовать анзац, где каждая функция устройства вычисляет одну конкретную реализацию ансамбля Лоренца.
struct lorenz_system { struct lorenz_functor { template< class T > __host__ __device__ void operator()( T t ) const { // unpack the parameter we want to vary and the Lorenz variables value_type R = thrust::get< 3 >( t ); value_type x = thrust::get< 0 >( t ); value_type y = thrust::get< 1 >( t ); value_type z = thrust::get< 2 >( t ); thrust::get< 4 >( t ) = sigma * ( y - x ); thrust::get< 5 >( t ) = R * x - y - x * z; thrust::get< 6 >( t ) = -b * z + x * y ; } }; lorenz_system( size_t N , const state_type &beta ) : m_N( N ) , m_beta( beta ) { } template< class State , class Deriv > void operator()( const State &x , Deriv &dxdt , value_type t ) const { thrust::for_each( thrust::make_zip_iterator( thrust::make_tuple( boost::begin( x ) , boost::begin( x ) + m_N , boost::begin( x ) + 2 * m_N , m_beta.begin() , boost::begin( dxdt ) , boost::begin( dxdt ) + m_N , boost::begin( dxdt ) + 2 * m_N ) ) , thrust::make_zip_iterator( thrust::make_tuple( boost::begin( x ) + m_N , boost::begin( x ) + 2 * m_N , boost::begin( x ) + 3 * m_N , m_beta.begin() , boost::begin( dxdt ) + m_N , boost::begin( dxdt ) + 2 * m_N , boost::begin( dxdt ) + 3 * m_N ) ) , lorenz_functor() ); } size_t m_N; const state_type &m_beta; };
<state_type><thrust::device_vector>или<device_vector>используется. Длина состояния3N, гдеN— число систем. Система кодируется в этот вектор таким образом, что всеxкомпоненты приходят первыми, затем каждыйукомпонентов и, наконец, каждыйzкомпонентов. Внедрение функции устройства является простой задачей, вам нужно только разложить кортеж, происходящий от итераторов зип.
Кроме системы без возмущений, нам необходимо вычислить систему, включающую линейные уравнения, управляющие временной эволюцией малых возмущений. Используя метод сверху, это просто, с небольшой сложностью, что кортежи Thrust имеют максимальную проходимость 10. Но это лишь небольшая проблема, так как мы можем создать итератор с итераторами. Таким образом, итератор зип верхнего уровня содержит один итератор зип для состояния, один нормальный итератор для параметра и один итератор зип для производного. Доступ к элементам этого кортежа в функции системы тогда прост, вы распаковываете кортеж<thrust::get<>()>. Мы не будем показывать код здесь, он большой. Его можно найти здесьи легко понять.
Кроме того, нам нужен наблюдатель, который определяет норму возмущений, нормализует их и усредняет логарифм нормы. Определяется функтор устройства, который используется в этом наблюдателе.
struct lyap_functor { template< class T > __host__ __device__ void operator()( T t ) const { value_type &dx = thrust::get< 0 >( t ); value_type &dy = thrust::get< 1 >( t ); value_type &dz = thrust::get< 2 >( t ); value_type norm = sqrt( dx * dx + dy * dy + dz * dz ); dx /= norm; dy /= norm; dz /= norm; thrust::get< 3 >( t ) += log( norm ); } };
Обратите внимание, что этот функтор манипулирует состоянием, то есть возмущениями.
Сейчас мы завершаем весь код для вычисления Ляпуновских экспонентов. Во-первых, мы должны определить вектор состояния. Этот вектор содержит6Nзаписи, состояниеx,y,zи его возмущенияdx,dy,dz. Мы начинаем их так, чтоx=y=z=10,dx=1иdy=dz=0. Определяем степперный тип, управляемый Runge-Kutta Dormand-Prince 5 степпер. Мы начинаем с некоторой интеграции, чтобы преодолеть переходное поведение. Для этого мы не вовлекаем возмущение и запускаем алгоритм только на состоянииx,y,zбез какого-либо наблюдателя. Обратите внимание, какBoost.Rangeиспользуется для частичной интеграции вектора состояния без возмущений (первая половина всего состояния). После переходного периода интегрирована полная система с возмущениями и рассчитаны и записаны Ляпуновские экспоненты<stdout>.
state_type x( 6 * N ); // initialize x,y,z thrust::fill( x.begin() , x.begin() + 3 * N , 10.0 ); // initial dx thrust::fill( x.begin() + 3 * N , x.begin() + 4 * N , 1.0 ); // initialize dy,dz thrust::fill( x.begin() + 4 * N , x.end() , 0.0 ); // create error stepper, can be used with make_controlled or make_dense_output typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type; lorenz_system lorenz( N , beta ); lorenz_perturbation_system lorenz_perturbation( N , beta ); lyap_observer obs( N , 1 ); // calculate transients integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz , std::make_pair( x.begin() , x.begin() + 3 * N ) , 0.0 , 10.0 , dt ); // calculate the Lyapunov exponents -- the main loop double t = 0.0; while( t < 10000.0 ) { integrate_adaptive( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , lorenz_perturbation , x , t , t + 1.0 , 0.1 ); t += 1.0; obs( x , t ); } vector< value_type > lyap( N ); obs.fill_lyap( lyap ); for( size_t i=0 ; i<N ; ++i ) cout << beta_host[i] << "\t" << lyap[i] << "\n";
Полный пример можно найти наlorenz_parameters.cu.
Статья Using CUDA (or OpenMP, TBB, ...) via Thrust раздела Chapter 1. Boost.Numeric.Odeint Tutorial может быть полезна для разработчиков на c++ и boost.
реклама |
