Разбиение на блоки в C++ AMP

Visual Studio 11 сделает мейнстримовой поддержку гетерогенных вычислений за счет технологии C++ Accelerated Massive Parallelism (C++ AMP). Я ознакомил вас с этой технологией в предыдущей статье в этом номере, которую я считаю подготовительной для чтения данной статьи. Поэтому, если вы еще не прочитали ее, пожалуйста, сделайте это сейчас.

Прежде чем описывать методику оптимизации при программировании GPU под названием «разбиение на блоки» (tiling), хочу напомнить, что в предыдущей статье вы узнали об index<N>, extent<N>, array_view<T,N>, restrict(amp) и parallel_for_each. И можете использовать C++ AMP API для реализации собственных параллельных алгоритмов обработки данных, например перемножения матриц, приведенного в предыдущей статье и повторяемого здесь на рис. 1.

Рис. 1. Перемножение матриц, простая модель

1  void MatMul(vector<int>& vC, const vector<int>& vA,
     const vector<int>& vB, int M, int N, int W )
2  {
3    array_view<const int, 2> a(M, W, vA), b(W, N, vB);
4    array_view<int, 2> c(M, N, vC);
5    c.discard_data();
6    parallel_for_each(c.extent, [=](index<2> idx) restrict(amp)
7    {
8      int row = idx[0]; int col = idx[1];
9      int sum = 0;
10     for(int i = 0; i < b.extent[0]; i++)
11       sum += a(row, i) * b(i, col);
12     c[idx] = sum;
13   });
14   c.synchronize();
15 }

Если вы не знакомы с перемножением матриц, почитайте материалы по ссылке bit.ly/HiuP.

В этой статье я расскажу о разбиении на блоки (tiling) — наиболее распространенном способе оптимизации при программировании графических процессоров (GPU). Единственная причина введения такой оптимизации в ваш алгоритм — дополнительный уровень производительности, которого вы можете добиться за счет повторного использования данных и более эффективных шаблонов доступа к памяти. Разбиение на блоки обеспечивает эффективное использование иерархии памяти GPU на более низком уровне, чем это возможно при использовании простой модели, известной вам из первой вводной статьи.

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

Класс tiled_extent

Вы знаете, что parallel_for_each принимает объект extent в качестве первого аргумента. Этот объект описывает домен вычислений (compute domain), т. е. сколько потоков (размер) и какой формы (число измерений) будут выполнять вычисления. Рассмотрим следующие два примера extent:

extent<1> e(12);   // 12 потоков в одномерном пространстве
extent<2> ee(2,6); // 12 потоков в двухмерном пространстве

Для parallel_for_each существует перегруженная версия с поддержкой разбиения на блоки, которая принимает аргумент tiled_extent. Этот объект описывает, как разбить исходный extent на блоки равного размера. Вы можете получить tiled_extent вплоть до ранга, равного трем. Если у вас больше измерений, придерживайтесь простой модели или проведите рефакторинг кода. Общее количество потоков в блоке (tile) не может превышать 1024.

Самый простой способ получить объект tiled_extent — вызвать непараметризованный шаблонный метод tile в extent, который возвращает tiled_extent для данного объекта extent. Чтобы получить соответствующие объекты tiled_extent для предыдущих двух примеров extent, можно написать нечто в таком духе:

tiled_extent<6> t_e = e.tile<6>();        // e.rank==t_e.rank
tiled_extent<2,2> t_ee = ee.tile<2, 2>(); // ee.rank==t_ee.rank

Соответствующую схему см. на рис. 2, где показано, как extent разбивается на меньшие подмножества данных, которые в терминологии C++ AMP называются блоками (tiles).

Рис. 2. Объект tiled_extent является объектом extent, разделенным на блоки

Числа, которые вы выбираете для передачи в качестве аргументов-шаблонов (template arguments), должны быть известны на этапе компиляции. Кроме того, они должны равномерно разделять измерения глобального extent, передаваемого в parallel_for_each:

• e[0]=12 можно разделить с помощью t_e.tile_extent[0]=6;
• ee[0]=2 можно разделить с помощью t_ee.tile_extent[0]=2;
• ee[1]=6 можно разделить с помощью t_ee.tile_extent[1]=2.

По соображениям производительности минимальный размер блока в наименьшем значащем измерении (least-significant dimension) должен быть равен 16. Варьирование размера блока может увеличивать или ухудшать производительность в зависимости от используемого вами оборудования. Мой совет: не пытайтесь делать этого — просто выбирайте такое число, которое обеспечивает одинаковую производительность на широком спектре оборудования, начиная с 16 и пробуя кратные этому значению четные числа.

Класс tiled_index

Вы знаете, что при вызове parallel_for_each передается лямбда с вашим кодом в качестве второго аргумента. Ваш код получает объект index, и вы можете считать этот объект идентификатором потока. Например:

parallel_for_each(my_extent, [=] (index<2> idx) restrict (amp) {
  my_array_view[idx] = ...
});

Когда вы разбиваете на блоки extent, передаваемый в parallel_for_each, сигнатура лямбды принимает tiled_index. Объект tiled_index состоит из четырех объектов index. Например, вы можете по-прежнему получить объект index, ожидаемый вами, через свойство объекта tiled_index:

parallel_for_each (my_extent.tile<2, 2>() , [=] (tiled_index<2, 2> t_idx) restrict (amp) {
  my_array_view [t_idx.global] = ...
});

При написании алгоритмов, использующих разбиение на блоки (блочных алгоритмов), вам, вероятно, потребуется знать локальный индекс в блоке (а не только глобальный индекс в домене вычислений в целом). Вы можете получить этот объект index через локальное свойство объекта tiled_index. В некоторых алгоритмах полезно знать индекс блока относительно других блоков, участвующих в вычислениях, а также глобальный индекс источника этого блока. Вы можете обращаться к соответствующим объектам index через свойства tile и tile_origin объекта tiled_index.

Используя двухмерный extent из предыдущего примера (extent<2>(2,6).tile<2,2>()), вы можете увидеть значения ранее упомянутых свойств tiled_index для выделенного на рис. 3 квадрата.

Рис. 3. Свойства tiled_index, возвращающие объекты index

Возвращаемся к перемножению матриц, но с разбиением на блоки

На рис. 1 вы видели реализацию перемножения матриц с использованием простой модели C++ AMP. Как изменить этот алгоритм, чтобы можно было явным образом применить разбиение на блоки с помощью средств, о которых вы уже знаете (используя tiled_extent и tiled_index)?

Решение показано на рис. 4 (изменения по сравнению с предыдущим листингом выделены полужирным).

Рис. 4. Перемножение матриц с разбиением на блоки без применения tile_static

3  array_view<const int, 2> a(M, W, vA), b(W, N, vB);
4  array_view<int, 2> c(M, N, vC);
5  c.discard_data();
6  parallel_for_each(c.extent.tile<16,16>(),
     [=](tiled_index<16,16> t_idx) restrict(amp)
7  {
8  int row = t_idx.global[0]; int col = t_idx.global[1];
9  int sum = 0;
10 for(int i = 0; i < b.extent[0]; i++)
11   sum += a(row, i) * b(i, col);
12 c[t_idx.global] = sum;
13 });
14 c.synchronize();

В строке 6 я вызвал метод tile объекта extent, выбрав размер блока (16×16 в этом примере), и изменил лямбду так, чтобы она принимала tiled_index с подходящими под размеры блока аргументами-шаблонами. В теле лямбды я заменил все вхождения idx на t_idx.global (строки 8 и 12). Это механистическое преобразование — первое, что вы должны делать для всех своих алгоритмов под C++ AMP, когда решаете ввести в них поддержку разбиения на блоки. Это первый — но не единственный — шаг на пути от простой модели к блочной.

Стоит отметить еще раз, что при таком изменении вы должны обеспечить, чтобы выбранный размер блоков был кратным измерениям глобального extent. В моем примере предполагается использование квадратных матриц, в которых каждое измерение можно без остатка разделить на 16. Кроме того, обычной практикой является запись размера блока в статическую const-переменную типа int или аргумент-шаблон.

В примере с простой моделью перемножения матрицы на рис. 1 система «за кулисами» сама осуществляет разбиение на блоки. Поэтому вам не приходится заботиться о требованиях к корректности деления. Но вот что простая модель не может сделать за вас (а значит, вы должны делать это сами), — это необходимы этап 2 разбиения на блоки: изменение алгоритма так, чтобы он использовал память tile_static и, как правило, один или более других объектов index. Прежде чем углубиться в эту тематику, давайте совершим краткий экскурс в иерархию памяти графических процессоров.

Краткий обзор иерархии памяти GPU

В сети много материалов по характеристикам графических процессоров. Материалы по-разному подаются в зависимости от конкретной компании — разработчика оборудования и даже семейства графических процессоров. В этом разделе я предлагаю вам составить высокоуровневую и приближенную картину, независимую от конкретной аппаратной платформы (каждый поставщик оборудования использует свою терминологию, и в каждой аппаратной платформе есть свои тонкости), — т. е. исключительно с точки зрения разработчика, применяющего C++ AMP.

В GPU есть глобальная память, в которой содержатся ваши данные array и array_view. С каждым потоком также сопоставляются регистры, куда обычно и помещаются ваши локальные переменные (если только драйвер вашей видеокарты не поправляет ваши просчеты — например, когда вы пытаетесь использовать слишком много регистров для оборудования, на котором выполняется ваш код, он будет перебрасывать их содержимое в глобальную память). Доступ к глобальной памяти осуществляется гораздо медленнее, чем к регистрам; например, для доступа к регистру достаточно одного такта GPU в сравнении с 1000 тактами для доступа к глобальной памяти.

Кроме того, каждому обрабатывающему элементу GPU отводится небольшой буфер памяти (скажем, 16–48 Кб в зависимости от оборудования). Обращение к ним осуществляется значительно быстрее, чем к глобальной памяти; к примеру, около 10 тактов. Эта область памяти, также называемая локальным хранилищем данных, является программируемым кешем. Кеши центрального процессора управляются автоматически и прозрачно для вас, поэтому любой выигрыш в производительности становится доступен вам безо всяких усилий с вашей стороны. С графическими процессорами дело обстоит прямо противоположным образом: вы должны самостоятельно управлять этим кешем GPU, копируя в него данные и извлекая их оттуда. В некоторых моделях программирования этот кеш называют общей памятью (shared memory), в других — локальной (local memory), а в третьих — общей памятью группы (group shared memory). В C++ AMP этот программируемый кеш называют памятью tile_static (подробности потом).

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

Часть глобальной памяти доступна всем потокам, и срок жизни данных в ней превышает таковой время вычислений с помощью parallel_for_each, поэтому последующий вызов parallel_for_each может оперировать с теми же данными. Значение регистра доступно лишь одному потоку, и его срок жизни совпадает с временем существования потока. Часть памяти tile_static является общей для подмножества всех потоков, которая в терминологии C++ AMP называется блоком потоков (tile of threads), и ее срок жизни равер таковому для блока потоков. Теперь вы наверняка начинаете понимать, зачем нужно разбиение вашей вычислительной операции на блоки: не зная, к какому потоку относится тот или иной блок, вы не сможете использовать эту быструю память tile_static.

Память tile_static можно рассматривать как «арендуемую» блоком потоков до окончания его выполнения, после чего эту память занимает другой блок. Поэтому, с логической точки зрения, потоки из разных блоков изолированы друг от друга в отношении памяти tile_static.

Использование нового класса хранилища tile_static

Для доступа к памяти tile_static вы используете класс хранилища tile_static. Это второе усовершенствование, которое C++ AMP вносит в язык C++ (первое — restrict, о котором я рассказывал в предыдущей статье).

Вы можете использовать этот новый класс хранилища только в функциях restrict(amp), только в том случае, если parallel_for_each разбивается на блоки, и только для переменных, объявленных внутри этого блока функции. Указатели и ссылки нельзя помечать как tile_static, а любые неявные конструкторы/деструкторы переменных tile_static не вызываются. Как правило (но не всегда), ваши переменные tile_static представляют собой массивы, и обычно они пропорциональны размеру блока, например:

tile_static float t[16][16]; // обращаемся к t для доступа
                                       //к статической памяти блока

Наиболее распространенный способ использования преимуществ памяти tile_static — выявить те области глобальной памяти, к которым ваш алгоритм может обращаться более одного раза. Затем скопировать эти области в память tile_static (лишь раз заплатив высокую цену за обращение к глобальной памяти) и изменить алгоритм так, чтобы он использовал копию данных в tile_static (к которой можно повторно обращаться очень быстро), не требуя многократного доступа к этим данным в глобальной памяти. Программируемый кеш имеет малый размер, поэтому в переменные tile_static нельзя скопировать все данные array и array_view. Блочный алгоритм (с разбиением на блоки) обычно сложнее, так как приходится обходить эту проблему, копируя из глобальной памяти только данные размером, соответствующим блоку.

Прежде чем обсуждать более реалистичный пример, иллюстрирующий сказанное в предыдущем абзаце, рассмотрим упрощенный, надуманный и содержащий ошибки пример, сфокусированный исключительно на использовании tile_static, где я пытаюсь суммировать все элементы матрицы:

1  static const int TS = 2;
2  array_view<int, 2> av(2, 6, my_vector);
3  parallel_for_each(av.extent.tile<TS,TS>(),
     [=](tiled_index<TS,TS> t_idx) restrict(amp)
4  {       
5    tile_static int t[TS][TS];   
6    t[t_idx.local[0]][t_idx.local[1]] = av[t_idx.global];
7
8    if (t_idx.local == index<2>(0,0)) {
9      t[0][0] = t[0][0] + t[0][1] + t[1][0] + t[1][1];             
10     av[t_idx.tile_origin] = t[0][0];
11   }
12 });
13 int sum = av(0,0) + av(0,2) + av(0,4); // три tile_origins

Этот код можно охарактеризовать одним словом: ужас. В строках 9 и 13 используется тот факт, что размер блока равен 2×2=4 и что общий размер — 2×6=12 (соответственно у нас три блока), тогда как настоящие реализации следует писать так, чтобы изменение размера блока или extent в целом не требовало изменения любого другого кода. Производительность этого кода также ужасна, потому что в нем применяется наивный алгоритм и его ветвление не поддаются распараллеливанию. В этом алгоритме отсутствует повторное использование данных, а значит применение памяти tile_static ничем не оправданно. Кроме того, в этом примере есть ошибка корректности, о которой я расскажу позже. Однако, как бы ни был ужасен этот код, он достаточно прост для полного новичка в области операций с хранилищем tile_static, чтобы разобраться в данном коде — и значение имеет только это.

В строке 6 каждый поток в каждом блоке копирует данные в память tile_static из глобальной памяти. В строке 9 только первый поток каждого блока будет суммировать результаты для этого блока в первую позицию в памяти tile_static блока. Затем в строке 10 тот же поток (в каждом блоке) будет сохранять результат обратно в глобальной памяти. По окончании вычислений на ускорителе три суммы от трех блоков будут складываться в переменной sum потоком, выполняемым на центральном процессоре.

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

Класс tile_barrier

Состояние гонок (race condition) в предыдущем примере возникает между строками 6 и 9. В строке 9 поток с локальным индексом (0,0) предполагает, что все значения уже хранятся в tile_static-переменной t, но это справедливо, только если остальные потоки в блоке выполнили строку 6. Поскольку потоки обычно выполняются независимо друг от друга, на это предположение полагаться нельзя.

Здесь нужно каким-то образом указать, что строка 9 не должна выполняться, пока все потоки не в блоке не выполнят строку 6. Для выражения этого ограничения используется класс tile_barrier и один из его методов wait. Вы не можете сами конструировать объект tile_barrier, но можете получить его через свойство barrier из tiled_index, переданного в вашу лямбду. Поэтому состояние гонок можно устранить, вставив после строки 6:

7 t_idx.barrier.wait();

Заметьте, что вы не смогли бы поместить эту строку до строки 9 и внутри условного блока. Вызов tile_barrier::wait должен появляться в том месте, где все потоки в блоке либо достигают этой точки, либо все вместе обходят ее. Если бы барьер разрешалось помещать в другие места, то в случае, когда один из потоков не выполнил wait, программа зависла бы, ожидая этот поток.

Многие из таких условий для гонок компилятор умеет сам помечать флагом, остальные поможет найти отладчик, если в Visual Studio 11 вы включите параметр GPU Memory Access Exceptions в диалоге Exceptions, выбрав из меню Debug команду Exceptions.

Пример: блочное перемножение матриц

Помните пример с перемножением матрицы на рис. 4? Теперь пора заняться второй частью упражнения с разбиением на блоки, где вы также задействуете преимущества памяти tile_static memory и будете использовать локальную индексацию и объект tile_barrier.

Прежде чем показать решение, напомню, что на моем лэптопе перемножение матриц с применением простой модели C++ AMP ускоряется более чем в 40 раз по сравнению с использованием обычного последовательного кода на центральном процессоре для случая M=N=W=1024. Решение с разбиением на блоки, которое вы вскоре увидите, с TS=16 (а значит, 16×16 = 256 потоков на блок) работает еще в два раза быстрее, чем на основе простой модели! Этот показатель учитывает передачу данных, поэтому, если вы уже передали данные и провели ряд вычислений с данными, ускорение выполнения должно быть гораздо больше, чем в два раза, по сравнению с версией, где не используется память tile_static. Насколько сложнее должен быть алгоритм, чтобы получить такой прирост производительности?

Полный код блочного перемножения матриц показан на рис. 5.

Рис. 5. Перемножение матриц с разбиением на блоки и применением tile_static

1  void MatMul(vector<int>& vC, const vector<int>& vA,
     const vector<int>& vB, int M, int N, int W )
2  {
3    array_view<const int,2> a(M, W, vA), b(W, N, vB);
4    array_view<int,2> c(M, N, vC);  
5    c.discard_data();
6
7    parallel_for_each(c.extent.tile<TS,TS>(),
       [=](tiled_index<TS,TS> t_idx) restrict(amp) 
8    {
9      int row = t_idx.local[0]; int col = t_idx.local[1];
10     tile_static int locA[TS][TS], locB[TS][TS];
11     int sum = 0;
12     for (int i = 0; i < a.extent[1]; i += TS) {
13       locA[row][col] = a(t_idx.global[0], col + i);
14       locB[row][col] = b(row + i, t_idx.global[1]);
15       t_idx.barrier.wait();
16
17       for (int k = 0; k < TS; k++)
18         sum += locA[row][k] * locB[k][col];           
19       t_idx.barrier.wait();
20     }
21     c[t_idx.global] = sum;
22   });
23   c.synchronize();
24 }

Хотя код на рис. 5 работает для блоков и экстентов любых размеров (при соблюдении описанных ранее правил), проще всего разобраться в нем, используя малые размеры. Такие размеры не обеспечивают должной производительности, но помогают ухватить суть происходящего. Поэтому допустим, что блоки имеют размер 2×2 (TS=2), а M=2, N=6 и W=4. Эта конфигурация приведена на рис. 6.

Рис. 6. Пример перемножения матриц с использованием 12 потоков в блоках 2×2

Чтобы вникнуть в код, следите только за одним блоком (вторым из трех) и только за одним потоком — тем, который вычисляет результат 160 (он выделен на рис. 6; выделенная область включает строку A и столбец B, к которым нужно обращаться этому потоку, чтобы подсчитывать результаты и помещать их в ячейку C).

Давайте пройдем код на рис. 5, поглядывая на полезную картинку на рис. 6.

Когда i=0, в окне Parallel Watch на рис. 7 показываются значения релевантных переменных вплоть до вложенного цикла в строке 16.

Рис. 7. Когда i=0, значения других переменных показываются в каждом столбце; при этом каждая строка является потоком

Как видите, четыре потока блока совместно копируют данные в два tile_static-массива locA и locB из матриц A и B. После того как все они встречаются на барьере в строке 15, они входят во вложенный цикл на строке 17. Значения соответствующих переменных при k=0, а затем при k=1 см. на рис. 8.

Рис. 8. По-прежнему i=0, значения других переменных показываются для k=0 и k=1

В этой точке поток, за которым вы следите, вычисляет частичную сумму 1×4 и на второй итерации переменной k добавляет ее к 2×10, что дает сумму, равную 24 (рис. 8). Далее он встречается с другими потоками на барьере в строке 19. Теперь они готовы войти во вторую итерацию внешнего цикла. Заметьте, что переменная i будет иметь значение 2, и на рис. 9 показываются новые значения переменных вплоть до этого барьера.

Рис. 9. Когда i=2, значения других переменных показываются в этой таблице

И вновь четыре потока блока совместно копируют данные в два tile_static-массива locA и locB из матриц A и B, двигаясь вправо по матрице A и вниз по матрице B. После того как все они встречаются на барьере в строке 15, они вновь входят во вложенный цикл на строке 17; Значения соответствующих переменных при k=0, а затем при k=1 см. на рис. 10.

В этой точке, согласно рис. 10, поток, за которым вы следите, добавляет к 24 результат произведения 3×16, а на второй итерации k прибавляет еще и результат умножения 4×22, что в итоге дает сумму, равную 160. Затем он встречается с остальными потоками на барьере в строке 19. Потоки готовы войти в третью итерацию внешнего цикла, но обнаруживают, что условие продолжения цикла (loop condition) больше не отвечает переменной i, поэтому они пропускают код до строки 21. Поток, за которым вы следите, использует свой глобальный индекс для обновления глобальной памяти в C.

Рис. 10. По-прежнему i=2, значения других переменных показываются для k=0 и k=1

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

Чтобы вам было легче понять этот код, я выбрал размер блока, равный 4 (TS=2), и малые экстенты, но в реальном коде они будут больше. Результаты по производительности, которыми я поделился с вами, достигаются при 16×16=256 потоков в каждом блоке, поэтому я повторно использовал порцию данных 256 раз из быстрой памяти вместо того, чтобы каждый раз обращаться к ним в глобальной памяти. Блочная реализация перемножения матриц дает дополнительный выигрыш от более эффективных шаблонов доступа к памяти для матрицы A (операции с матрицами B и C выполняются эффективно во всех реализациях), но смешивание памяти (memory coalescing) выходит за рамки этой статьи.

Вот и все, что я хотел рассказать о блочном перемножении матриц, использующим преимущества памяти tile_static. Если вы хотите побольше попрактиковаться в работе с блочными алгоритмами, см. примеры в блоге группы Parallel Programming in Native Code по ссылке bit.ly/xZt05R.

Выбор оптимального варианта

В этой статье я познакомил вас с разбиением на блоки — наиболее распространенным способом оптимизации кода, использующего C++ AMP.

Вы научились использовать три новых класса C++ AMP (tiled_extent, tiled_index и tile_barrier), а также работать с классом хранилища tile_static. Теперь вы знаете, что можно начать с реализации на основе простой модели (без явного разбиения на блоки). Потом можно модифицировать реализацию, вызывая функцию tile объекта extent (для выбора размера блока) и изменив сигнатуру лямбды так, чтобы она принимала tiled_index (с аргументами-шаблонами для того же размера блока).

На этом чисто механистический этап 1 завершается. Для второго (последнего) этапа вы должны переписать свой алгоритм под использование памяти tile_static с соответствующей синхронизацией через tile_barrier. Этот этап требует творческого подхода, где для каждого алгоритма приходится создавать совершенно новое решение. Простая и блочная реализации перемножения матриц демонстрируют уровень сложности, вводимый разбиением на блоки, а также показывают, почему их применение может давать колоссальное увеличение производительности. Выбор оптимального варианта остается за вами.

Стоит упомянуть, что такое программирование кеша может давать выигрыш и в коде, выполняемом на центральном процессоре, поскольку это поможет более эффективной работе системе автоматического прозрачного управления кешем. Кроме того, ваш код, разбитый на блоки, станет более подходящим для векторизации (распараллеливания) интеллектуальными компиляторами.

Помимо двух моих статей, в блоге группы Parallel Programming (bit.ly/bWfC5Y) есть множество ресурсов по C++ AMP (документация, советы, шаблоны, ссылки на видеоролики и примеры), которые я настоятельно рекомендую изучить.