Анализ влияния размера рабочей группы на производительность OpenCL-реализации вычислительного алгоритма на примере метода Гаусса решения СЛАУ

Информационные технологии 236 Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 5 (2), с. 236-241

УДК 519.684.4

АНАЛИЗ ВЛИЯНИЯ РАЗМЕРА РАБОЧЕЙ ГРУППЫ НА ПРОИЗВОДИТЕЛЬНОСТЬ OPENCL-РЕАЛИЗАЦИИ ВЫЧИСЛИТЕЛЬНОГО АЛГОРИТМА НА ПРИМЕРЕ МЕТОДА ГАУССА РЕШЕНИЯ СЛАУ

© 2012 г. Д.А. Баранов, И.В. Влацкая

Оренбургский госуниверситет

mois@mail.osu.ru

Поступила в редакцию 10.09.2012

Описываются некоторые особенности стандарта ОрепСЬ, проявляющиеся при использовании на аппаратном обеспечении архитектуры СИОА. В частности, анализируется влияние размера локальной рабочей группы на время вычислений. В качестве примера используется ОрепСЬ -реализация метода Гаусса решения СЛАУ.

Ключевые слова: GPGPU, ОрепСЬ, SIMT, параллельные вычисления, метод Гаусса, СЛАУ.

Современную информационную инфраструктуру сложно представить без многопроцессорных вычислительных комплексов. Параллельные вычисления используются для удовлетворения большинства информационных потребностей человека, начиная от теоретических исследований до сугубо практических приложений, таких как прогнозирование погоды. Не удивительно, что растущая потребность в эффективной обработке больших объёмов данных привела к появлению большого количества подходов, технологий и, как следствие, программно-аппаратных архитектур для исполнения параллельных алгоритмов обработки данных. Особый интерес представляет относительно новое направление «General-Purpose Graphics Processing Units» (GPGPU) - программирование графического процессора общего назначения. GPGPU появилось в результате активного развития графических процессоров в последнее десятилетие под влиянием игровой и кино индустрии. На данный момент графические процессоры обладают значительно большей теоретической пиковой производительностью, чем центральные процессоры той же ценовой категории.

До появления GPGPU решение задач общего назначения на графических процессорах было крайне затруднительно. Объясняется это особенностью архитектуры современных графических процессоров, оптимизированной для решения задач обработки геометрических и графических данных. Ситуация изменилась с появлением разработок AMD (на тот момент ATI) и NVIDIA для программирования графических процессоров собственного производства: ATI

Stream SDK и NVIDIA CUDA (Compute Unified Device Architecture). Несмотря на то, что графические процессоры AMD и NVIDIA решали одинаковые задачи, их технологии GPGPU несовместимы и работают исключительно на аппаратном обеспечении фирмы-производителя. В последствии группа Khronos, в которую вошли AMD, NVIDIA и другие производители графических процессоров и вычислительных устройств, разработала стандарт GPGPU OpenCL, целью которого была поддержка большого количества платформ. Следует отметить, что AMD позже отказалась от собственной разработки в пользу OpenCL. На данный момент OpenCL работает на графических процессорах как AMD, так и NVIDIA, а также на любых x86 и amd64-совместимых центральных процессорах. Поэтому OpenCL является универсальным инструментом разработки высокопроизводительных параллельных алгоритмов, обеспечивающим работу большей части современной компьютерной техники.

Тем не менее, стандарт OpenCL относительно молод (первая версия была опубликована 9 декабря 2008), и на данный момент существуют плохо документированные особенности, которые отрицательно влияют на его применимость. Так, анализ источников информации по данной тематике показал, что многие программисты сталкиваются с проблемами производительности OpenCL, особенно на аппаратном обеспечении производства NVIDIA. Существуют результаты экспериментов, подтверждающие эти проблемы, и результаты сравнительных испытаний с другими технологиями GPGPU, в частности CUDA. Причём некоторые авторы явно

заостряют внимание на преимуществе технологии СИБЛ над ОреиСЬ в производительности. В качестве примера на рис. 1 приведён график зависимости времени расчётов клеточных автоматов игры «Жизнь» от размеров ребра кубического поля из работы [1].

30

25

20

15

10

5

OpenCL x-xCUDA

4 32 64

128

256

Рис. 1. Зависимость времени расчетов от размера ребра кубического «поля» при использовании CUDA и OpenCL

Кроме того, некоторые результаты экспериментов показывают наличие «выбросов» времени вычислений в зависимости от объёма входных данных [2].

Для подробного изучения причины появления «выбросов» времени вычислений необходимо рассмотреть важные особенности архитектуры графических процессоров NVIDIA и стандарта OpenCL. Следует заметить, что как бы не называлась аппаратная архитектура того или иного поколения графических процессоров NVIDIA (например Fermi, Kepler, Tesla), фактически все они являются подмножеством архитектуры CUDA.

CUDA является проприетарной (закрытой) программно-аппаратной архитектурой и, следовательно, пригодна лишь для графических процессоров NVIDIA. В то же время, OpenCL является открытым стандартом, разработанным группой Khronos, в состав которой входит множество производителей вычислительных устройств. Стандарт OpenCL разработан для использования на большом количество платформ и устройств, в т.ч. на CPU и GPU. Следовательно, область его применения значительно шире. Тем не менее, архитектура CUDA оказала большое влияние на разрабатываемый стандарт, поскольку к тому времени CUDA уже заняла большую часть рынка высокопроизводительных вычислений.

Основой архитектуры CUDA является масштабируемый массив многопоточных мультипроцессоров (Streaming Multiprocessors, SMs). При запуске вычислений CUDA разделяет вычислительные потоки на блоки и распределяет их по мультипроцессорам. На одном мульти-

процессоре может выполняться более одного блока процессов, но один блок процессов не может одновременно выполняться на нескольких мультипроцессорах. Мультипроцессоры разработаны для одновременного выполнения большого количество потоков, для этого используется специальная архитектура, называемая SIMT (Single-Instruction, Multiple-Thread -одна команда, множество потоков). Кроме того, используется конвейер команд для достижения параллелизма на уровне команд в пределах одного потока. Однако следует помнить, что ядра графического процессора устроены гораздо проще CPU, в частности, отсутствует предсказание ветвлений. Более подробно архитектура CUDA описана в [3, 4].

Для проведения вычислительного эксперимента, необходимого для анализа проблемы выбросов, использовался метода Гаусса [5] решения СЛАУ. Рассмотрим простейшую реализацию этого метода на OpenCL.

Реализация алгоритма метода Гаусса с использованием OpenCL разделена на четыре функции («ядра» в терминологии OpenCL), по два на прямой и обратный ход. На каждом шаге прямого хода вызывается сначала первое (gauss_fwd_pre), затем второе ядро (gauss_fwd). Пусть текущая итерация имеет номер к, размерность СЛАУ - п. В первом ядре прямого хода проводится вычисление делителей для каждой строки на текущей итерации:

4+1 =A-,i=k+\, n.

Ak

Akk

Во втором ядре прямого хода вычисляются новые значения матрицы А и вектора Ь:

4+1 = A - 4Ak, i=k+1,n,j=k+1,n

bk+1 = bkAk, i = k + 1, n.

Заметим, что формула (2) не обнуляет элементы к-го столбца, это необходимо, так как все

лк+1

элементы Л^ вычисляются параллельно и

используют значения к-го столбца. Обнуление этого столбца должно производиться после вы, к+1

числения всех элементов Л , однако это не

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

Рассмотрим обратный ход метода Гаусса. На каждой итерации обратного хода вызываются два ядра: gauss_bwd_prepare и gauss_bwd. В первом ядре обратного хода вычисляется хк:

x-bL Xk Ak.

Akk

Очевидно, для первого шага эта формула работает, а для последующих это становится возможно, благодаря второму ядру:

Ь*+1 = Ь* - хк4к, 1= й-1.

Как видно из (4), вычисленное значение хк

сразу используется для модификации следующих уравнений. Фактически, происходит перенос компонентов хкЛк в правую часть уравнений, но в левой части они не обнуляются, а просто игнорируются в дальнейшем. Ядро (3), очевидно, неоптимально, поскольку использует всего один процессор, но его использование более оправдано, чем пересылка из памяти графического процессора в оперативную память и обратно.

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

Описываемые ниже наблюдения были получены при использовании первого подхода, так как второй подход значительно усложняет анализ экспериментальных данных. Следует отметить, что при использовании данной реализации алгоритма отмечаются временные выбросы в виде локальных максимумов на матрицах некоторой размерности, т. е. значительное падение производительности. График зависимости времени выполнения алгоритма от размерности матрицы изображён на рис. 2.

Подробное исследование проблемы показало, что «выбросы» встречаются на матрицах, для которых выполняется условие:

3 !к,[[ + 1,к] = 0,к ф 1,

где N - размерность СЛАУ; [а, Ь] - остаток от деления а на Ь; N + 1 - размерность СЛАУ + вектор-столбец В.

t, c

3000 -

2000

1000

0

2000

4000

6000

8000

10000

Рис. 2. Зависимость времени выполнения алгоритма от размерности матрицы

Фактически, выражение (5) означает, что N + 1 - простое число. Анализ рекомендаций по программированию для CUDA от NVIDIA [3] позволил сделать предположение о существовании проблемы оптимального разделения глобальной рабочей области на рабочие группы (следует заметить, что в описываемой реализации размер рабочий группы не задавался явно и, следовательно, устанавливался драйвером). Действительно, согласно документации, размер рабочей группы должен быть делителем размера глобальной рабочей области и при этом не превышать некоторого ограничения конкретного вычислительного устройства (величина ограничения может быть получена с помощью функции OpenCL clGetDevicelnfo, вызванной с параметром CL_DEVICE_MAX_WORK_GROUP _SIZE). Для большинства графических процессоров это ограничение составляет 1024 - это максимальное число вычислительных процессов в составе одной группы, оно составляется из произведения чисел всех рабочих процессов по всем измерениям внутри группы. Кроме того, особенности архитектуры графических процессоров, в частности концепция «Warp» от NVIDIA, позволяют максимизировать производительность параллельного алгоритма путём задания размеров рабочей группы кратными специальной величине, которая так же может быть запрошена у вычислительного устройства. Эта величина называется «Preferred Work Group Size Multiplier» (PWGSM) и её можно получить с помощью функции OpenCL clGetKernelWork GroupInfo, вызванной с параметром CL_KER-NEL_PREFERRED_WORK_GROUP_SIZE_MU-LTIPLE. Для графических процессоров NVIDIA эта величина равна 32. Фактически, «Warp» -это минимальная группа вычислительных процессов (до 32), которой может управлять планировщик графического процессора.

Практика подтвердила, что в случае простых чисел в качестве размерности глобальной рабочей области, драйвер автоматически устанавливает размер рабочей группы, равным единст-

венному делителю размера глобальной рабочей области, меньшему 1024, т.е. 1. Очевидно, что большое число рабочих групп, состоящих из 1 элемента, не могут быть эффективно распределены по вычислительным устройствам графического процессора. В этом случае большая часть времени работы алгоритма уходит на работу планировщика и сводит на нет преимущества Б1МТ-архитектуры.

Заметим, что многие авторы уделяют большое внимание работе с памятью графических процессоров, например [6]. Однако обычно опускаются важные вопросы разделения глобальной рабочей области на рабочие группы.

Для решения проблемы выбросов было проведено исследование производительности описанной выше реализации метода Гаусса в зависимости от размера рабочей группы и размерности СЛАУ.

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

[SZ1,OGS• К] = 0, К = OSG},

SZ2 = М - OSG • К,

где М - размерность матрицы (М = N + 1 в описанном примере); OSG - оптимальный размер группы; SZ1 - размер первой части глобальной рабочей области; SZ2 - размер второй части глобальной рабочей области; {а, Ь} - целая часть от деления а на Ь.

Как видно, размер второй части является остатком от деления и может быть произвольным числом меньше OSG. Способ вычисления второй части, в общем случае, зависит от специфики задачи и может производиться как отдельной реализацией алгоритма (например последовательной СРИ-версией), так и повторным вызовом OpenCL-ядра с установкой соответствующего размера глобальной рабочей области (для этой цели можно создавать отдельную копию ядра). Очевидно, величина остатка SZ2 будет влиять на общее время работы алгоритма, поэтому оптимальное значение OSG потенциально может не привести к наиболее быстрому достижению результата.

Ниже приведены результаты экспериментов

с различными размерами локальной рабочей группы, устанавливаемыми для второго ядра описанной выше реализации. Размер матрицы так же варьируется.

На рис. 3 приведён график зависимости времени вычислений от размера рабочей группы для размерностей СЛАУ от 100 до 900 с шагом 100.

мс 16000

12000

8000

4000

0

Рис. 3. Зависимость времени вычислений от размера рабочей группы

Как видно из графика 3, размер группы, равный единице, максимально неэффективен на любом размере матрицы. Рассмотрим данный график более детально, для этого необходимо отбросить результат с размером группы, равным единице. Результат приведён на рис. 4.

Рис. 4. Зависимость времени вычислений от размера рабочей группы, секция 32-896

Из графика 4 следует, что для большинства размерностей матрицы оптимальный размер группы лежит в пределах 192-384 элемента. Тем не менее, на матрицах размерности < 1000 остаток SZ2 может оказывать значительное влияние на результат, поэтому необходимо рассмотреть матрицы большей размерности. На рис. 5 изображён график зависимости времени вычислений от размера рабочей группы для размерностей матрицы от 1000 до 10000 с шагом 1000.

Рассмотрим секцию 128-384 графика 5 подробнее, она приведена на рис. 6.

Согласно графику 6, оптимальным размером рабочей группы является 192, что подтверждается табличными данными (с подробными ре-

N = 100 N = 200 N = 300 -N = 400 N = 500 N = 600 N = 700 -N = 800 N = 900

1 64

192

320 448

576

704 832

зультатами экспериментов, а также с исходным кодом можно ознакомиться в открытом репози-тории на ОкЬиЬ [7]).

Рис. 5. Зависимость времени вычислений от размера рабочей группы для матриц большой размерности

Рис. 6. Зависимость времени вычислений от размера рабочей группы для матриц большой размерности, секция 128-384

Очевидно, значение 192 не является оптимальным для любого вычислительного устройства. Тем не менее, это значение достаточно просто выводится из характеристик устройства, которые можно получить через функции OpenCL. Характеристики вычислительного устройства (в данном случае графического процессора), на котором проводился данный эксперимент, можно получить на сайте производителя [8]. Однако там не указан один важный параметр - число мультипроцессоров (Comput Unit в терминологии OpenCL). Как уже упоминалось, один вычислительный блок (рабочая группа) не может быть разделён на несколько мультипроцессоров, но несколько блоков могут выполняться на одном мультипроцессоре. Число Compute Unit вычислительного устройства можно получить с помощью вызова функции clGetDevicelnfo с параметром CL_DEVICE_ MAX_COMPUTE_UNITS.

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

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

_ „ totalCoreCount

mpCoreCount =-,

mpCount

где totalCoreCount - общее число ядер («потоки») устройства, в данном случае 144; mpCount - число мультипроцессоров, для данного вычислительного устройства это число составляет 3.

Следовательно, для данного устройства число ядер мультипроцессора составляет 144/3 = = 48. Наименьшее общее кратное 32 и 48, большее 144 равно 192 - это и есть оптимальный размер рабочей группы для данного вычислительного устройства, что подтверждается результатами эксперимента.

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

Результаты данной работы используются для разработки OpenCL-реализации параллельного алгоритма обучения нечёткой нейронной сети Такаги - Сугено - Канга.

Работа выполнена при поддержке федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. (контракт 14.B37.21.0176).

Список литературы

1. Алексеенко А.Е., Казённов А.М. Реализация клеточных автоматов «игра «Жизнь» с применением технологий CUDA и OpenCL // Математические основы и численные методы моделирования. 2010. Т. 2. № 3. С. 323-326. МФТИ.

2. Запрягаев С.А., Карпушин А.А. Применение графического процессора в задаче многократного численного интегрирования на основе обобщённой формулы Симпсона // Вестник Воронежского госуниверситета. 2011. №1. С. 150-156.

3. NVIDIA CUDA C programming Guide, 2012.

4. Запрягаев С. А., Карпушин А. А. Применение графического процессора в ресурсоемких вычис-

лениях на базе библиотеки OpenCL // Вестник Воронежского госуниверситета. 2010. №2.

5. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.

6. Манушин И.А. Использование технологии OpenCL для разработки высоконагруженных при-

ложений // RSDN Magazine. 2012. № 1. С. 12-20. К-Пресс.

7. https://github.com/Dem0n3D/solecl

8. http://www.nvidia.ru/object/product-geforce-gt-555 m-ru.html

ANALYSIS OF THE EFFECT OF WORK-GROUP SIZE ON OPENCL ALGORITHM PERFORMANCE AS EXEMPLIFIED BY SLAE SOLUTION BY THE GAUSS METHOD

D.A. Baranov, I.V. Vlatskaya

This paper describes some features of the OpenCL standard that are manifested when using the CUDA hardware architecture. In particular, the effect is analyzed of the size of the local work- group on the computation time. As an example, we give here an OpenCL implementation of the SLAE solution by the Gauss method.

Keywords: GPGPU, OpenCL, SIMT, parallel computing, Gauss method, SLAE.