Использование GPU для вычислений с помощью C++ AMP
125C# и .NET Framework --- Оптимизация приложений .NET Framework --- Использование GPU для вычислений с помощью C++ AMP
До сих пор в обсуждении приемов параллельного программирования мы рассматривали только ядра процессора. Мы приобрели некоторые навыки распараллеливания программ по нескольким процессорам, синхронизации доступа к совместно используемым ресурсам и использования высокоскоростных примитивов синхронизации без применения блокировок.
Однако, существует еще один способ распараллеливания программ - графические процессоры (GPU), обладающие большим числом ядер, чем даже высокопроизводительные процессоры. Ядра графических процессоров прекрасно подходят для реализации параллельных алгоритмов обработки данных, а большое их количество с лихвой окупает неудобства выполнения программ на них. В этой статье мы познакомимся с одним из способов выполнения программ на графическом процессоре, с использованием комплекта расширений языка C++ под названием C++ AMP.
Расширения C++ AMP основаны на языке C++ и именно поэтому в данной статье будут демонстрироваться примеры на языке C++. Однако, при умеренном использовании механизма взаимодействий в. NET, вы сможете использовать алгоритмы C++ AMP в своих программах для .NET. Но об этом мы поговорим в конце статьи.
Введение в C++ AMP
По сути, графический процессор является таким же процессором, как любые другие, но с особым набором инструкций, большим количеством ядер и своим протоколом доступа к памяти. Однако между современными графическими и обычными процессорами существуют большие отличия, и их понимание является залогом создания программ, эффективно использующих вычислительные мощности графического процессора.
Современные графические процессоры обладают очень маленьким набором инструкций. Это подразумевает некоторые ограничения: отсутствие возможности вызова функций, ограниченный набор поддерживаемых типов данных, отсутствие библиотечных функций и другие. Некоторые операции, такие как условные переходы, могут стоить значительно дороже, чем аналогичные операции, выполняемые на обычных процессорах. Очевидно, что перенос больших объемов кода с процессора на графический процессор при таких условиях требует значительных усилий.
Количество ядер в среднем графическом процессор значительно больше, чем в среднем обычном процессоре. Однако некоторые задачи оказываются слишком маленькими или не позволяют разбивать себя на достаточно большое количество частей, чтобы можно было извлечь выгоду от применения графического процессора.
Поддержка синхронизации между ядрами графического процессора, выполняющими одну задачу, весьма скудна, и полностью отсутствует между ядрами графического процессора, выполняющими разные задачи. Это обстоятельство требует синхронизации графического процессора с обычным процессором.
Сразу возникает вопрос, какие задачи подходят для решения на графическом процессоре? Имейте в виду, что не всякий алгоритм подходит для выполнения на графическом процессоре. Например, графические процессоры не имеют доступа к устройствам ввода/вывода, поэтому у вас не получится повысить производительность программы, извлекающей ленты RSS из интернета, за счет использования графического процессора. Однако на графический процессор можно перенести многие вычислительные алгоритмы и обеспечить массовое их распараллеливание. Ниже приводится несколько примеров таких алгоритмов (этот список далеко не полон):
увеличение и уменьшение резкости изображений, и другие преобразования;
быстрое преобразование Фурье;
транспонирование и умножение матриц;
сортировка чисел;
инверсия хеша «в лоб».
Отличным источником дополнительных примеров может служить блог Microsoft Native Concurrency, где приводятся фрагменты кода и пояснения к ним для различных алгоритмов, реализованных на C++ AMP.
C++ AMP - это фреймворк, входящий в состав Visual Studio 2012, дающий разработчикам на C++ простой способ выполнения вычислений на графическом процессоре и требующий лишь наличия драйвера DirectX 11. Корпорация Microsoft выпустила C++ AMP как открытую спецификацию, которую может реализовать любой производитель компиляторов.
Фреймворк C++ AMP позволяет выполнять код на графических ускорителях (accelerators), являющихся вычислительными устройствами. С помощью драйвера DirectX 11 фреймворк C++ AMP динамически обнаруживает все ускорители. В состав C++ AMP входят также программный эмулятор ускорителя и эмулятор на базе обычного процессора, WARP, которые служит запасным вариантом в системах без графического процессора или с графическим процессором, но в отсутствие драйвера DirectX 11, и использует несколько ядер и инструкции SIMD.
А теперь приступим к исследованию алгоритма, который легко можно распараллелить для выполнения на графическом процессоре. Реализация ниже принимает два вектора одинаковой длины и вычисляет поточечный результат. Сложно представить что-либо более прямолинейное:
void VectorAddExpPointwise(float* first, float* second, float* result, int length)
{
for (int i = 0; i < length; ++i) {
result[i] = first[i] + exp(second[i]);
}
}Чтобы распараллелить этот алгоритм на обычном процессоре, требуется разбить диапазон итераций на несколько поддиапазонов и запустить по одному потоку выполнения для каждого из них. Мы посвятили достаточно много времени в предыдущих статьях именно такому способу распараллеливания нашего первого примера поиска простых чисел - мы видели, как можно это сделать, создавая потоки вручную, передавая задания пулу потоков и используя Parallel.For и PLINQ для автоматического распараллеливания. Вспомните также, что при распараллеливании похожих алгоритмов на обычном процессоре мы особо заботились, чтобы не раздробить задачу на слишком мелкие задания.
Для графического процессора эти предупреждения не нужны. Графические процессоры имеют множество ядер, выполняющих потоки очень быстро, а стоимость переключения контекста значительно ниже, чем в обычных процессорах. Ниже приводится фрагмент, пытающийся использовать функцию parallel_for_each из фреймворка C++ AMP:
#include <amp.h>
#include <amp_math.h>
using namespace concurrency;
void VectorAddExpPointwise(float* first, float* second, float* result, int length)
{
array_view <const float,1> avFirst (length, first);
array_view <const float,1> avSecond(length, second);
array_view <float,1> avResult(length, result);
avResult.discard_data();
parallel_for_each(avResult.extent, [=](index<1> i) restrict(amp)
{
avResult[i] = avFirst[i] + fast_math::exp(avSecond[i]);
});
avResult.synchronize();
}Теперь исследуем каждую часть кода отдельно. Сразу заметим, что общая форма главного цикла сохранилась, но первоначально использовавшийся цикл for был заменен вызовом функции parallel_for_each. В действительности, принцип преобразования цикла в вызов функции или метода для нас не нов - ранее уже демонстрировался такой прием с применением методов Parallel.For() и Parallel.ForEach() из библиотеки TPL.
Далее, входные данные (параметры first, second и result) обертываются экземплярами array_view. Класс array_view служит для обертывания данных, передаваемых графическому процессору (ускорителю). Его шаблонный параметр определяет тип данных и их размерность. Чтобы выполнить на графическом процессоре инструкции, обращающиеся к данным, первоначально обрабатываемым на обычном процессоре, кто-то или что-то должен позаботиться о копировании данных в графический процессор, потому что большинство современных графических карт являются отдельными устройствами с собственной памятью. Эту задачу решают экземпляры array_view - они обеспечивают копирование данных по требованию и только когда они действительно необходимы.
Когда графический процессор выполнит задание, данные копируются обратно. Создавая экземпляры array_view с аргументом типа const, мы гарантируем, что first и second будут скопированы в память графического процессора, но не будут копироваться обратно. Аналогично, вызывая discard_data(), мы исключаем копирование result из памяти обычного процессора в память ускорителя, но эти данные будут копироваться в обратном направлении.
Функция parallel_for_each принимает объект extent, определяющий форму обрабатываемых данных и функцию для применения к каждому элементу в объекте extent. В примере выше мы использовали лямбда-функцию, поддержка которых появилась в стандарте ISO C++2011 (C++11). Ключевое слово restrict (amp) поручает компилятору проверить возможность выполнения тела функции на графическом процессоре и отключает большую часть синтаксиса C++, который не может быть скомпилирован в инструкции графического процессора.
Параметр лямбда-функции, index<1> объекта, представляет одномерный индекс. Он должен соответствовать используемому объекту extent - если бы мы объявили объект extent двумерным (например, определив форму исходных данных в виде двумерной матрицы), индекс также должен был бы быть двумерным. Пример такой ситуации приводится чуть ниже.
Наконец, вызов метода synchronize() в конце метода VectorAddExpPointwise гарантирует копирование результатов вычислений из array_view avResult, произведенных графическим процессором, обратно в массив result.
На этом мы заканчиваем наше первое знакомство с миром C++ AMP, и теперь мы готовы к более подробным исследованиям, а так же к более интересным примерам, демонстрирующим выгоды от использования параллельных вычислений на графическом процессоре. Сложение векторов - не самый удачный алгоритм и не самый лучший кандидат для демонстрации использования графического процессора из-за больших накладных расходов на копирование данных. В следующем подразделе будут показаны два более интересных примера.
Умножение матриц
Первый «настоящий» пример, который мы рассмотрим, - умножение матриц. Для реализации мы возьмем простой кубический алгоритм умножения матриц, а не алгоритм Штрассена, имеющий время выполнения, близкое к кубическому ~O(n2.807). Для двух матриц: матрицы A размером m x w и матрицы B размером w x n, следующая программа выполнит их умножение и вернет результат - матрицу C размером m x n:
void MatrixMultiply(int* A, int m, int w, int* B, int n, int* C)
{
for (int i = 0; i < m; ++i)
{
for (int j = 0; j < n; ++j)
{
int sum = 0;
for (int k = 0; k < w; ++k)
{
sum += A[i * w + k] * B[k * w + j];
}
C[i*n + j] = sum;
}
}
}Распараллелить эту реализацию можно несколькими способами, и при желании распараллелить этот код для выполнения на обычном процессоре правильным выбором был бы прием распараллеливания внешнего цикла. Однако графический процессор имеет достаточно большое количество ядер и распараллелив только внешний цикл, мы не сможем создать достаточное количество заданий, чтобы загрузить работой все ядра. Поэтому имеет смысл распараллелить два внешних цикла, оставив внутренний цикл нетронутым:
void MatrixMultiply (int* A, int m, int w, int* B, int n, int* C)
{
array_view <const int,2> avA(m, w, A);
array_view <const int,2> avB(w, n, B);
array_view <int,2> avC(m, n, C);
avC.discard_data();
parallel_for_each (avC.extent, [=](index <2> idx) restrict(amp)
{
int sum = 0;
for (int k = 0; k < w; ++k)
{
sum + = avA(idx[0]*w, k) * avB(k*w, idx[1]);
}
avC[idx] = sum;
});
}Эта реализация все еще близко напоминает последовательную реализацию умножения матриц и пример сложения векторов, приводившиеся выше, за исключением индекса, который теперь является двумерным и доступен во внутреннем цикле с применением оператора []. Насколько эта версия быстрее последовательной альтернативы, выполняемой на обычном процессоре? Умножение двух матриц (целых чисел) размером 1024 х 1024 последовательная версия на обычном процессоре выполняет в среднем 7350 миллисекунд, тогда как версия для графического процессора - держитесь крепче - 50 миллисекунд, в 147 раз быстрее!
Моделирование движения частиц
Примеры решения задач на графическом процессоре, представленные выше, имеют очень простую реализацию внутреннего цикла. Понятно, что так будет не всегда. В блоге Native Concurrency, ссылка на который уже приводилась выше, демонстрируется пример моделирования гравитационных взаимодействий между частицами. Моделирование включает бесконечное количество шагов; на каждом шаге вычисляются новые значения элементов вектора ускорений для каждой частицы и затем определяются их новые координаты. Здесь распараллеливанию подвергается вектор частиц - при достаточно большом количестве частиц (от нескольких тысяч и выше) можно создать достаточно большое количество заданий, чтобы загрузить работой все ядра графического процессора.
Основу алгоритма составляет реализация определения результата взаимодействий между двумя частицами, как показано ниже, которую легко можно перенести на графический процессор:
// здесь float4 - это векторы с четырьмя элементами,
// представляющие частицы, участвующие в операциях
void bodybody_interaction (
float4& acceleration, const float4 p1, const float4 p2) restrict(amp)
{
float4 dist = p2 – p1;
// w здесь не используется
float absDist = dist.x*dist.x + dist.y*dist.y + dist.z*dist.z;
float invDist = 1.0f / sqrt(absDist);
float invDistCube = invDist*invDist*invDist;
acceleration + = dist*PARTICLE_MASS*invDistCube;
}Исходными данными на каждом шаге моделирования является массив с координатами и скоростями движения частиц, а в результате вычислений создается новый массив с координатами и скоростями частиц:
struct particle
{
float4 position, velocity;
// реализации конструктора, конструктора копирования и
// оператора = с restrict(amp) опущены для экономии места
};
void simulation_step(array <particle,1> & previous,
array <particle,1> & next, int bodies)
{
extent <1> ext(bodies);
parallel_for_each (ext, [&](index <1> idx) restrict(amp)
{
particle p = previous[idx];
float4 acceleration(0, 0, 0, 0);
for (int body = 0; body < bodies; ++body)
{
bodybody_interaction (acceleration, p.position, previous[body].position);
}
p.velocity + = acceleration*DELTA_TIME;
p.position + = p.velocity*DELTA_TIME;
next[idx] = p;
});
}С привлечением соответствующего графического интерфейса, моделирование может оказаться очень интересным. Полный пример, представленный командой разработчиков C++ AMP, можно найти в блоге Native Concurrency. На моей системе с процессором Intel Core i7 и видеокартой Geforce GT 740M, моделирование движения 10 000 частиц выполняется со скоростью ~2.5 кадра в секунду (шагов в секунду) с использованием последовательной версии, выполняющейся на обычном процессоре, и 160 кадров в секунду с использованием оптимизированной версии, выполняющейся на графическом процессоре - огромное увеличение производительности.
Прежде чем завершить этот раздел, необходимо рассказать еще об одной важной особенности фреймворка C++ AMP, которая может еще больше повысить производительность кода, выполняемого на графическом процессоре. Графические процессоры поддерживают программируемый кеш данных (часто называемый разделяемой памятью (shared memory)). Значения, хранящиеся в этом кеше, совместно используются всеми потоками выполнения в одной мозаике (tile). Благодаря мозаичной организации памяти, программы на основе фреймворка C++ AMP могут читать данные из памяти графической карты в разделяемую память мозаики и затем обращаться к ним из нескольких потоков выполнения без повторного извлечения этих данных из памяти графической карты. Доступ к разделяемой памяти мозаики выполняется примерно в 10 раз быстрее, чем к памяти графической карты. Иными словами, у вас есть причины продолжить чтение.
Чтобы обеспечить выполнение мозаичной версии параллельного цикла, методу parallel_for_each передается домен tiled_extent, который делит многомерный объект extent на многомерные фрагменты мозаики, и лямбда-параметр tiled_index, определяющий глобальный и локальный идентификатор потока внутри мозаики. Например, матрицу 16x16 можно разделить на фрагменты мозаики размером 2x2 (как показано на рисунке ниже) и затем передать функции parallel_for_each:
extent <2> matrix(16,16);
tiled_extent <2,2> tiledMatrix = matrix.tile <2,2> ();
parallel_for_each (tiledMatrix, [=](tiled_index <2,2> idx)
restrict(amp)
{
// ...
});
Каждый из четырех потоков выполнения, принадлежащих одной и той же мозаике, могут совместно использовать данные, хранящиеся в блоке.
При выполнении операций с матрицами, в ядре графического процессора, взамен стандартного индекса index<2>, как в примерах выше, можно использовать idx.global. Грамотное использование локальной мозаичной памяти и локальных индексов может обеспечить существенный прирост производительности. Чтобы объявить мозаичную память, разделяемую всеми потоками выполнения в одной мозаике, локальные переменные можно объявить со спецификатором tile_static.
На практике часто используется прием объявления разделяемой памяти и инициализации отдельных ее блоков в разных потоках выполнения:
parallel_for_each(tiledMatrix, [=](tiled_index <2,2> idx) restrict(amp)
{
// 32 байта совместно используются всеми потоками в блоке
tile_static int local[2][2];
// присвоить значение элементу для этого потока выполнения
local[idx.local[0]][idx.local[1]] = 42;
});Очевидно, что какие-либо выгоды от использования разделяемой памяти можно получить только в случае синхронизации доступа к этой памяти; то есть, потоки не должны обращаться к памяти, пока она не будет инициализирована одним из них. Синхронизация потоков в мозаике выполняется с помощью объектов tile_barrier (напоминающего класс Barrier из библиотеки TPL) - они смогут продолжить выполнение только после вызова метода tile_barrier.Wait(), который вернет управление только когда все потоки вызовут tile_barrier.Wait. Например:
parallel_for_each (tiledMatrix, [](tiled_index <2,2> idx) restrict(amp)
{
// 32 байта совместно используются всеми потоками в блоке
tile_static int local[2][2];
// присвоить значение элементу для этого потока выполнения
local[idx.local[0]][idx.local[1]] = 42;
// idx.barrier - экземпляр tile_barrier
idx.barrier.wait();
// Теперь этот поток может обращаться к массиву "local",
// используя индексы других потоков выполнения!
});
Теперь самое время воплотить полученные знания в конкретный пример. Вернемся к реализации умножения матриц, выполненной без применения мозаичной организации памяти, и добавим в него описываемую оптимизацию. Допустим, что размер матрицы кратен числу 256 - это позволит нам работать с блоками 16 х 16. Природа матриц допускает возможность поблочного их умножения, и мы можем воспользоваться этой особенностью (фактически, деление матриц на блоки является типичной оптимизацией алгоритма умножения матриц, обеспечивающей более эффективное использование кеша процессора).
Суть этого приема сводится к следующему. Чтобы найти Ci,j (элемент в строке i и в столбце j в матрице результата), нужно вычислить скалярное произведение между Ai,* (i-я строка первой матрицы) и B*,j (j-й столбец во второй матрице). Однако, это эквивалентно вычислению частичных скалярных произведений строки и столбца с последующим суммированием результатов. Мы можем использовать это обстоятельство для преобразования алгоритма умножения матриц в мозаичную версию:
void MatrixMultiply(int* A, int m, int w, int* B, int n, int* C)
{
array_view <const int,2> avA(m, w, A);
array_view <const int,2> avB(w, n, B);
array_view <int,2> avC(m, n, C);
avC.discard_data();
parallel_for_each (avC.extent.tile <16,16> (), [=](tiled_index <16,16> idx) restrict(amp)
{
int sum = 0;
int localRow = idx.local[0], localCol = idx.local[1];
for (int k = 0; k <w; k += 16)
{
tile_static int localA[16][16], localB[16][16];
localA[localRow][localCol] = avA(idx.global[0], localCol + k);
localB[localRow][localCol] = avB(localRow + k, idx.global[1]);
idx.barrier.wait();
for (int t = 0; t <16; ++t)
{
sum + = localA[localRow][t]*localB[t][localCol];
}
// чтобы избежать затирания разделяемой памяти в следующей итерации
idx.barrier.wait();
}
avC[idx.global] = sum;
});
}Суть описываемой оптимизации в том, что каждый поток в мозаике (для блока 16 х 16 создается 256 потоков) инициализирует свой элемент в 16 х 16 локальных копиях фрагментов исходных матриц A и B. Каждому потоку в мозаике требуется только одна строка и один столбец из этих блоков, но все потоки вместе будут обращаться к каждой строке и к каждому столбцу по 16 раз. Такой подход существенно снижает количество обращений к основной памяти.
Чтобы вычислить элемент (i,j) в матрице результата, алгоритму требуется полная i-я строка первой матрицы и j-й столбец второй матрицы. Когда потоки мозаике 16x16, представленные на диаграмме и k=0, заштрихованные области в первой и второй матрицах будут прочитаны в разделяемую память. Поток выполнения, вычисляющий элемент (i,j) в матрице результата, вычислит частичное скалярное произведение первых k элементов из i-й строки и j-го столбца исходных матриц.
В данном примере применение мозаичной организации обеспечивает огромный прирост производительности. Мозаичная версия умножения матриц выполняется намного быстрее простой версии и занимает примерно 17 миллисекунд (для тех же исходных матриц размером 1024 х 1024), что в 430 быстрее версии, выполняемой на обычном процессоре!
Прежде чем закончить обсуждение фреймворка C++ AMP, нам хотелось бы упомянуть инструменты (в Visual Studio), имеющиеся в распоряжении разработчиков. Visual Studio 2012 предлагает отладчик для графического процессора (GPU), позволяющий устанавливать контрольные точки, исследовать стек вызовов, читать и изменять значения локальных переменных (некоторые ускорители поддерживают отладку для GPU непосредственно; для других Visual Studio использует программный симулятор), и профилировщик, дающий возможность оценивать выгоды, получаемые приложением от распараллеливания операций с применением графического процессора. За дополнительной информацией о возможностях отладки в Visual Studio обращайтесь к статье «Пошаговое руководство. Отладка приложения C++ AMP» на сайте MSDN.
Альтернативы вычислений на графическом процессоре в .NET
До сих пор в этой статье демонстрировались примеры только на языке C++, тем не менее, есть несколько способов использовать мощь графического процессора в управляемых приложениях. Один из способов - использовать инструменты взаимодействий, позволяющие переложить работу с ядрами графического процессора на низкоуровневые компоненты C++. Это решение отлично подходит для тех, кто желает использовать фреймворк C++ AMP или имеет возможность использовать уже готовые компоненты C++ AMP в управляемых приложениях.
Другой способ - использовать библиотеку, непосредственно работающую с графическим процессором из управляемого кода. В настоящее время существует несколько таких библиотек. Например, GPU.NET и CUDAfy.NET (обе являются коммерческими предложениями). Ниже приводится пример из репозитория GPU.NET GitHub, демонстрирующий реализацию скалярного произведения двух векторов:
[Kernel]
public static void MultiplyAddGpu(double[] a, double[] b, double[] c)
{
int ThreadId = BlockDimension.X * BlockIndex.X + ThreadIndex.X;
int TotalThreads = BlockDimension.X * GridDimension.X;
for (int ElementIdx = ThreadId; ElementIdx <a.Length; ElementIdx += TotalThreads)
{
c[ElementIdx] = a[ElementIdx] * b[ElementIdx];
}
}Я придерживаюсь мнения, что гораздо проще и эффективнее освоить расширение языка (на основе C++ AMP), чем пытаться организовывать взаимодействия на уровне библиотек или вносить существенные изменения в язык IL.
Итак, после того как мы рассмотрели возможности параллельного программирования в .NET и использованием GPU наверняка ни у кого не осталось сомнений, что организация параллельных вычислений является важным способом повышения производительности. Во многих серверах и рабочих станциях по всему миру остаются неиспользуемыми бесценные вычислительные мощности обычных и графических процессоров, потому что приложения просто не задействуют их.
Библиотека Task Parallel Library дает нам уникальную возможность включить в работу все имеющиеся ядра центрального процессора, хотя при этом и придется решать некоторые интереснейшие проблемы синхронизации, чрезмерного дробления задач и неравного распределения работы между потоками выполнения.
Фреймворк C++ AMP и другие многоцелевые библиотеки организации параллельных вычислений на графическом процессоре с успехом можно использовать для распараллеливания вычислений между сотнями ядер графического процессора. Наконец, имеется, неисследованная ранее, возможность получить прирост производительности от применения облачных технологий распределенных вычислений, превратившихся в последнее время в одно из основных направлений развития информационных технологий.
