CC BY
38
Информатика, вычислительная техника и управление УДК 519.67 DOI: 10.14529/cmse160307
ОЦЕНКА ЛОКАЛЬНОСТИ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ, РЕАЛИЗУЕМЫХ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
© 2016 г. Н.А. Лиходед, М.А. Полещук
Белорусский государственный университет (220030 Республика Беларусь, Минск, пр. Независимости, д. 4) E-mail: likhoded@bsu.by, poleschuma@bsu.by Поступила в редакцию: 02.03.2016
Исследуется задача получения блоков операций и потоков операций параллельного алгоритма, приводящих к меньшему числу обращений к глобальной памяти и к эффективному использованию параллельными потоками вычислений кэшей и разделяемой памяти графического процессора. Сформулированы и доказаны утверждения, позволяющие оценить объем коммуникационных операций, порождаемых альтернативными вариантами задания размеров блоков вычислений, а также минимизировать число промахов кэша за счет использования временной и пространственной локальности данных с учетом размера и длины строк кэша. Исследования конструктивны и допускают программную реализацию для практического использования.
Ключевые слова: параллельные вычисления, графический процессор, минимизация объема коммуникационных операций, временная локальность, пространственная локальность.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Лиходед Н.А., Полещук М.А. Оценка локальности параллельных алгоритмов, реализуемых на графических процессорах // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2016. Т. 5, № 3. С. 96-111. DOI: 10.14529/cmse160307.
Введение
Время решения задачи на современном компьютере во многом определяется локальностью реализуемого алгоритма — степенью использования памяти с быстрым доступом. В этой работе в качестве компьютера, на котором требуется реализовать параллельную версию алгоритма, будем рассматривать графические процессоры (GPU) — мощные многоядерные вычислительные устройства. При вычислениях на GPU быстрым является процесс обращения к разделяемой памяти мультипроцессора и к кэшам, но не обращение к глобальной памяти GPU.
Графический процессор выполняет множество параллельных потоков вычислений. Потоки объединяются в блоки, а блоки объединяются в сетку. Синхронизация между потоками разных блоков не предусмотрена. Каждый блок потоков выполняется атомарно на одном из мультипроцессоров графического процессора. Должны быть указаны блоки, которые могут выполняться мультипроцессорами одновременно и независимо друг от друга.
Целью этой работы является разработка метода, позволяющего получать блоки вычислений с меньшим числом обращений к глобальной памяти. Используется анализ информационных зависимостей, порождающих коммуникационные операции. Сформулиро-
Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии — 2016».
ваны и доказаны утверждения, позволяющие ранжировать параметры размера блоков на основе асимптотических оценок объема коммуникационных операций. В работе [1] подобные оценки получены с использованием более сложного (но и более приспособленного для автоматизации) математического аппарата. Известен способ получения более точной оценки количества элементов, к которым осуществляется доступ при выполнении операций блока вычислений, но практическое использование известного результата довольно трудоемкое и требует привлечения специализированных средств автоматизации [2, 3].
В работе также указаны условия наличия временной локальности и условия наличия пространственной локальности данных. Существование локальности позволяет эффективно использовать параллельными потоками вычислений кэши и разделяемую память. Исследование проводится аналитическими методами. Другой способ исследования локальности данных — построение профиля работы приложения с памятью [4].
Статья организована следующим образом. В разделе 1 приведены необходимые для дальнейшего изложения сведения о рассматриваемом классе алгоритмов, информационных зависимостях между операциями, получении блоков макроопераций. В разделе 2 оценивается объем коммуникационных операций чтения и операций записи, порождаемых разбиением итераций алгоритма на блоки вычислений, приводится пример установления соотношения размеров блоков для минимизации объема коммуникационных операций. В разделе 3 даются рекомендации по выбору размеров потоков вычислений для достижения наименьшего числа кэш-промахов и организации совместного использования потоками данных. В разделе 4 обсуждаются результаты вычислительных экспериментов. В заключении суммируются основные результаты, намечаются направления дальнейших исследований.
1. Предварительные сведения
Будем считать, что алгоритм задан последовательной программой линейного класса1 [5]. Основную вычислительную часть такой программы составляют циклические конструкции; границы изменения параметров циклов задаются неоднородными формами, линейными по совокупности параметров циклов и внешних переменных. Предполагается, что в гнезде циклов имеется К выполняемых операторов 5р и используется L массивов а1 размерностей VI. Область изменения параметров циклов (область итераций) для оператора 5р и размерность этой области обозначим соответственно Ур и Пр.
Вхождением (1,р,с) будем называть с-е вхождение массива а1 в оператор 5р. Индексы
элементов 1-го массива, связанных с вхождением (!,р,с), выражаются функцией Fl вида
(3) = ¥1 ^ + Gl + /' ^,
3Ц,...,^) е^, N е ^, ^ е 7*G1 ^ е -,е Г1,
где N е 7^ — вектор внешних переменных алгоритма, з — число внешних переменных.
1 В литературе используются также термины «аффинные гнезда циклов», «алгоритмы с аффинными зависимостями».
Выполнение оператора при конкретных значениях р и вектора параметров циклов 1 будем называть операцией и обозначать 5^(7). Пара вхождений (1,а,1) и (/,р,д) порождает истинную зависимость 5а(1)—^5^(1), если: 5а(1) выполняется раньше 5р(1); 5а(1) переопределяет (изменяет) элемент массива а, а £р(1) использует на вхождении (/,р,д) в качестве аргумента тот же элемент массива; между операциями £а(1) и 5^(1) этот элемент не переопределяется.
Зависимости (информационные связи) £а(1)—5р(1) между операциями будем характеризовать векторами (Г'в=1-1. Если размерности векторов 1 и I не совпадают, то 1-1 определяется как разность векторов меньшей из размерностей. Векторы часто называют векторами зависимостей: позволяет для операции 5р(1) найти операцию £а(1),
от которой 5р(1) зависит.
Пусть в гнезде циклов имеется 0 наборов выполняемых операторов. Под набором операторов будем понимать один или несколько операторов, окруженных одним и тем же множеством циклов. Операторы и наборы операторов линейно упорядочены расположением их в записи алгоритма. Обозначим: V8, 1 <8<0, — области изменения параметров циклов, окружающих наборы операторов, п8 — размерность области V8, число циклов, окружающих 8-й набор операторов. Будем считать, что если оператор 5р принадлежит набору операторов с номером 8е, то п8 =Яр.
Выделим блоки вычислений путем разбиения циклов, окружающих операторы.
Пусть т8 = / , м8 = тах 7 , 1 <С<п8 , — предельные значения изменения
С ^(л,л.с с ^{м^8
параметров циклов. Зададим размеры блоков вычислений натуральными числами г^,...,. Параметр г^ обозначает число значений параметра приходящихся на один
блок 8-го набора операторов. Число частей 0^, на которые при формировании блоков разбивается область значений параметра цикла, окружающего 8-й набор операторов, находится согласно 0? = [(М^8 - т^ +1)/ г8| ([•] обозначает ближайшее «сверху» целое число). Нумеровать блоки вычислений будем по каждой координате в пределах от 0 до
08 -1, 1 <с< п8.
Формализованным способом разбить множество операций алгоритма на блоки и потоки вычислений можно путем применения тайлинга — преобразования алгоритма для получения макроопераций [6, 7].
Пример 1. Рассмотрим алгоритм Флойда—Уоршелла поиска кратчайших путей между всеми парами вершин графа (рис. 1).
^ к = 1, п ^ 1 = 1, п ^ 3 = 1, п
а(1,3) = т1п (а ( i,j) , а(1,к) + а(к3))
enddo enddo
Рис. 1. Основная часть алгоритма Флойда—Уоршелла
В гнезде циклов имеется один выполняемый оператор ^ (один набор операторов) и
используется один массив а размерности 2. Область изменения параметров циклов (область итераций) V1=V1 Z3| 1<г<п, 1<т<п} для оператора S1 имеет размерность 3. Для матриц Fцiq на вхождениях (1$,д) имеем:
F1,1,1 F1,1,2
0 1 0 0 0 1
F = 1,1,3
0 1 0 1 0 0
F1,1,4 =
1 0 0 0 0 1
В рассматриваемом примере векторы зависимостей d а'в будем для наглядности помечать не номерами операторов, а элементами массивов, фигурирующими на порождающих зависимости вхождениях. Например, вектор da(1']>'а(1'к> порождается зависимостью между данными а(,') на вхождении (1,а,1) =(1,1,1) в левой части оператора Sl, и данными а^,^) на вхождении (/,р,д)=(1,1,3) в правой части оператора. Укажем итерации, порождающие зависимости, и векторы зависимостей:
• S1(k-1,i;j)^S1(k;i;j): данное а(^), вычисленное на итерации является аргументом а(i,j) для вычислений на итерации (к^^); da(,J■>'a(г=(1,0,0);
• S1(k-1,i,k)^■S1(k,i,j): а^,к), вычисленное на итерации (к-1,^к), является аргументом для вычислений на итерации (к^^); da(1']>'а(1'к> =(1,0^-к);
• S1(k-1,k;j)^S1(k;i;j): а(к^), вычисленное на итерации (к-1,к^), является аргументом для вычислений на итерации (k,i,j); da(JJ>'a=(1^^,^).
Отметим, что для фиксированного k все операции не зависят друг от друга. Кроме того, непосредственно из записи алгоритма следует, что данные а(,') не обновляются, если i=k или j=k.
Выделим блоки вычислений. Напомним, блоки вычислений должны выполняться атомарно, независимо друг от друга. С учетом зависимостей алгоритма, такие блоки проще всего задать, если положить г=1, г2 и Г3 — параметры. Всего будет ^х^хОз
двумерных (2D) блоков, где Q=n, Q= ющий вид:
, Q3
Блочный алгоритм имеет следу-
n
n
Г2
r
3
do к = 1, n
do ibl = 0, Q2 - 1
do jbl = 0, Q3 - 1
do i = 1 + ibl*r2, min ( (ibl + 1)*r2, n)
do j = 1 + jbl*r3, min(( jbl + 1)*r3, n)
a(i,j) = min(a(i,j), a(i,k) + a(k,j))
enddo
enddo
enddo
enddo
enddo
Рис. 2. Основная часть алгоритма Флойда—Уоршелла с выделенными 2D блоками 2. Оценка объема коммуникационных операций 2.1. Выражения для оценки объема операций чтения и записи
Обозначим M8 = max (-M® - m®) +1 — наибольшее число итераций циклов, участвующих в получении блоков. Для простоты записи будем использовать обозначение M без
индекса 8 и подразумевать набор операторов 8е при упоминании оператора 5р. Отметим, что величина 08г^ имеет порядок М, поэтому г8, 08 могут быть величинами порядка М.
Определим термин «фиксированное данное массива» как конкретное, неизмененное содержимое соответствующей ячейки памяти. Введем в рассмотрение величины р^ ,
Рй$д, которые определим следующим образом: пусть число фиксированных данных, ис-
РП
пользуемых на вхождении (/,р,д), оценивается величиной О(М 1'р'д), а число фиксированных данных, используемых на вхождении (/,р,д) при фиксированном значении цикла
Рй £ -1
с параметром у оценивается величиной О(МК1') . Наличие индекса d в обозначениях подчеркивает, что если вхождение (/,р,д) порождает истинную зависимость 5а(2)^5р(1), то р1 ,р ,д и зависят от вектора Найти р^р л и р^д не представляет труда, если
детально понимать реализуемый алгоритм. Формализованный способ получения величин Р^р , р^д изложен в работе [1]. Этот способ использует функции Fl,р ,д, а также функции, описывающие информационную структуру алгоритма — покрывающие функции графа алгоритма [5] (называемые еще ^преобразованиями [8]).
Отметим, что число фиксированных данных, используемых (в пространстве размерности Яр-1) на вхождении при фиксированном значении цикла с параметром по порядку либо равно (р^д -1=Р^р д ), либо на 1 меньше (р^д -1=р^ д -1) числа фиксированных данных, используемых (в пространстве размерности Лр) при всех значениях цикла с параметром
Обозначим через й^'™* максимальные по модулю значения компонент векторов Сделаем естественное для практики предположение: компоненты векторов
по
модулю или не превосходят нескольких единиц, или неограниченно возрастают с ростом размера задачи. В первом случае ^(а'р'тах также не превосходят нескольких единиц, во
втором случае являются большими, значительно превосходящими г^8 .
Теорема 1. Пусть вхождение (/,р,?) порождает истинную зависимость: определение некоторого данного происходит на вхождении (1,а,1) в левой части оператора 5а, а использование — на вхождении (/,р,д) в правой части оператора 5р, причем в окружении операторов 5а и 5р имеется цикл с параметром Тогда, при получении блоков вычислений для реализации алгоритма на графическом процессоре, разбиение итераций цикла jz порождает коммуникационные операции чтения и записи, объем которых имеет следующие оценки:
• О08МР1 вл 1) операций чтения и операций записи, если 0 < аа'в'тах < г8 •
С '
ОР1 -р-д) операций чтения и О(МР1 р>д) операций записи, если ^а'Р'тах > Г
Р1 , в , д ? Р1, в , д •
ГМ Л Лр<а,в,д\ " ^а.В.тах ^ 8 ,С
О(М ) операций чтения и операций записи, если а^ > г^ , р1 р = р1 р ;
• не требуется ни операций чтения, ни операций записи, если d(a'в'max =0.
В случае, когда вхождение (1,р,с) не порождает истинной зависимости (происходит обращение к входным данным) или цикл с параметром имеется в окружении только
оператора £р, оценки следующие:
• О(О?Мргр-4 > операций чтения, если р^д Ф Р^ д;
Р d С
О(М г д д > операций чтения, если р1 £ = р
d
I ,р, д = р1 ,р, д '
Доказательство. Утверждения теоремы оценивают объем коммуникационных операций, порождаемых разбиением итераций цикла в зависимости от значений £-й координаты вектора и величин р^д, р/р д. Напомним, множество итераций цикла с параметром jz разбивается на О? частей, каждая часть имеет некоторый номер , 0 < < О? — 1, и содержит г^ итераций (кроме, возможно, части с номером О? — 1).
Пусть 0 < dIa'в'max < г?, определение данного происходит при выполнении операции (I(ц,... ,/Иа>> , а использование — при выполнении операции ^,••• ,У^» .
Число фиксированных данных, используемых на вхождении (1,р,с) при фиксированном значении цикла с параметром (т.е. в пространстве размерности Яр-1) оценивается
Рd—1
величиной О(М г 'р'д > . Тогда число фиксированных данных, используемых на вхождении (1,р,д) при всех значениях цикла с параметром (в пространстве размерности пв)
в
оценивается величиной —М—О(М р^д—1> = О(М р^д > (напомним, величина dZI'в'max не пре-
da'P'max ^ ' ^ ' ^
вышает нескольких единиц). С другой стороны, оценка числа таких данных есть
р^ d С
О(М ''р'д > . Таким образом, в рассматриваемом случае верно ргРд = р
d
1,в,д = р1,в,д .
В каждой из О? частей итераций цикла число фиксированных данных, которые
■г!
определяются при одном значении ' , но используются при другом, оценивается велика Дтах^/л гр/вд—1 \ /аДтах
чиной dz О(М > , причем независимо от значения г^ величина d^ , как было
условлено, не превышает нескольких единиц. Суммарный объем коммуникационных операций чтения и операций записи по всем частям определяется оценкой
О?О(М—1> = О (о?М —1>.
тт ^аАтах ^ ? ^ „ d,С / „ ^ / „ d,С 1 „ d \
Пусть dz > г, . Если р1 р дФр1 р (т.е. р1 рд -1=р1 р ), то число используемых на вхождении (1$^) фиксированных данных в каждой части и во всех О? частях оценива-
рd
ется одинаковой величиной О(М 1 'р'? >. Этой величиной следует оценить объем коммуникационных операций записи. Оценка объема коммуникационных операций чтения по
р1 а рd
мых фиксированных данных в каждой из О? частей оценивается величиной
всем О? частям есть О?О(Мрг,м > = О(О?Мрг"м > . Если = рdвд , то число используе-
<
О(г?МРг 1> . Суммарный объем коммуникационных операций чтения и операций запи-
еа
у частям определяется оценкой
б?0 (г?М рё^-1) = О (б?г?М рём-1) = 0(М ).
Если ё аДтах =0, то любое фиксированное данное определяется и используется при одном и том же значении параметра цикла . Коммуникационных операций не требуется.
Пусть вхождение (/,р,?) не порождает истинной зависимости (происходит обращение к входным данным) или порождает истинную зависимость, но цикл с параметром %
имеется в окружении только оператора 5р. Если р^д Ф рёр д (т.е. р^д -1=рё в д ), то объем коммуникационных операций (только чтение) по всем б? частям оценивается вели-
чиной б?0(Мр',в д) = О(б?Мр'в д) . Если р^д = рё в д , то объем коммуникационных операций (только чтение) по всем б? частям определяется оценкой
б?0 (г?М р<ё*л-1) = 0(М р'вд)
□
Доказательство теоремы для случая, когда информационная структура алгоритма представлена покрывающими функциями графа алгоритма, приведено в работе [1].
Пример 1 (продолжение). Напомним векторы зависимостей, порождаемые вхождениями массива а в правую часть оператора, и укажем оценки числа фиксированных данных, используемых на вхождениях.
Для второго вхождения а(,') (использование прежнего значения обновляемого элемента) ёа (г' 7а (г' 7) =(1,0,0); р11 2 =3 — число фиксированных данных, используемых на
вхождении, оценивается величиной О(п3); р^ = р^ = р^ =3 — число фиксированных данных, используемых на вхождении при фиксированном значении одного из циклов с параметрами %2=% Зз=к, оценивается величинами 0(МР1,1,2 *) =О(п).
Для третьего вхождения а^,к) (использование при фиксированном к на итерациях
элементов одного столбца) ёа(г' 7-1'а('к-1 =(1,0,%-к), рё 1 3 =2, р^ =2, рё 123 =2, рё 133 =3. Для четвертого вхождения а(к,%) (использование при фиксированном к на итерациях (^,к) элементов одной строки) ёа(г'7-1'а(к'7-1 =(1,-к,0), р1ё1,4 =2, рё14 =2, рё^ =3, р^ =2.
Оценим, используя теорему 1, объем коммуникационных операций, порождаемых разбиением итераций циклов.
Пусть С=1. Тогда для всех вхождений ё^'3'тах = ё^'3 =1. Имеет место случай
0 < ё^'3'тах < Г?, требуется О(б?Мр'-р-д 1) операций чтения и операций записи: для второго вхождения О^-уП ) операций, для третьего и четвертого вхождений — О^п) операций.
Пусть С=2. Для второго и третьего вхождений ё^' втах = ё™'3 =0, поэтому не требуется ни операций чтения, ни операций записи. Для четвертого вхождения (случай ё"' в'тах > г? , рЦр^дФрЦр ,д, (^=(1,1,4)) требуется О(б?МРёвд) =О(^2п2) операций чте-
рё 2 ния и 0(М^д) =О( п ) операций записи.
Пусть С=3. Для второго и четвертого вхождений d(a'в'max = d(a'в =0, не требуется ни операций чтения, ни операций записи. Для третьего вхождения (случай d(a'в'mаx > Г?, / А Ф ^
р?,в^дФр^^рд, (*,М=(1,1,3)) требуется О(0„п2) операций чтения и О(п) операций записи.
2.2. Ранжирование размеров блоков вычислений
Утверждения теоремы 1 позволяют найти асимптотику суммарного объема коммуникационных операций, порождаемых разбиением множества итераций ' на О? частей, и, следовательно, дают возможность ранжировать параметры размера блоков вычислений для минимизации объема коммуникационных операций. При ранжировании (установлении соотношений размеров относительно друг друга) параметров размера блоков вычислений параллельного алгоритма, реализуемого на графическом процессоре, следует в оценках объема коммуникационных операций принимать во внимание только оценки, содержащие множитель ; следует также учитывать, что величины г^ и
связаны обратно пропорциональной зависимостью.
Пример 1 (продолжение). Для ранжирования параметров Гр г2 и Г3 размера блоков вычислений, выполняемых на графическом процессоре, укажем, следуя результатам предыдущего подраздела, оценки объема коммуникационных операций, содержащие множитель О? — число частей, на которые при формировании блоков разбивается область значений параметра
Если С=1, то требуется 0(0^) операций чтения и операций записи для второго вхождения и О(0^п) операций чтения и операций записи для третьего и четвертого вхождений. Если С=2, то требуется О(02п ) операций чтения для четвертого вхождения. Если С=3, то требуется 0(0311) операций чтения для третьего вхождения.
Таким образом, сравнивая оценки коммуникационных издержек и по чтению, и по записи, приходим к выводу, что параметр размера блока вычислений г1 увеличивать
более выгодно, чем параметры т2 и Т3. В то же время, наиболее простой блочный алгоритм Флойда—Уоршелла (рис. 2) содержит 2D блоки, для которых ^=1. Блочный алгоритм с 3D блоками (рис. 3), для которых т1=т2=т3=т>1 получен в работе [9]. Представляет также интерес разработка блочного алгоритма с 3D блоками, для которых
ТХ>Т2, Т1>Тд.
do к = 1 + kbl*r, min ( (kbl + 1)*r, n)
do i = 1 + ibl*r, min ( (ibl + 1)*r, n)
do j = 1 + jbl*r, min ( (jbl + 1)*r, n)
a ( i,j) = min(a(i,j), a(i,k) + a(k,j)) enddo enddo enddo
Рис. 3. 3D блоки блочного алгоритма Флойда—Уоршелла
3. Условия наличия локальности данных
Потоки одного блока выполняются на мультипроцессоре частями-пулами, называемыми варп. Варп содержит до 32 потоков (два полуварпа по 16 потоков). Если несколько потоков одного полуварпа обращаются к одному и тому же банку разделяемой памяти, константной памяти или кэша текстурной памяти для доступа к различным данным, то происходит конфликт. Но конфликта не происходит, если несколько потоков одного полуварпа обращаются к одному и тому же данному (broadcast). В этом случае запрашиваемое данное передается только один раз, поэтому трафик сокращается в 16 раз (для архитектуры Fermi в 32 раза, так как объединение запросов в память происходит для 32 потоков).
Текстурная память может быть использована эффективно, если каждый поток по-луварпа запрашивает близко расположенные в памяти данные, то есть если доступ к памяти характеризуется пространственной локальностью. Наличие пространственной локальности позволяет эффективно использовать кэш текстурной памяти и разделяемую память параллельными потоками.
Зададим в блоках вычислений потоки вычислений посредством выделения блоков второго уровня. Размеры блоков второго уровня зададим натуральными числами ^ 2, r22, ..., rn,2. Параметр r^ обозначает число значений параметра j^, приходящихся на один поток (на один блок второго уровня). Обозначим через Q52 число частей, на которые при формировании потоков разбивается область значений цикла с параметром j^, приходящихся на один блок; Q^2=\Г'% /Г»21 • Выбор размеров r^2, оптимальных по числу
кэш-промахов, определяется в значительной степени наличием временной локальности и пространственной локальности данных.
Введем обозначения:
F/p — столбец с номером 5 матрицы F^
— номер ненулевого столбца матрицы F^q, правее которого находятся только нулевые столбцы матрицы; если самый правый столбец не нулевой, то по определению 5 =Пр;
— множество номеров линейно независимых столбцов матрицы F;„q (если таких множеств более одного, то зафиксировать любое);
F 'р q — матрица, составленная из столбцов матрицы F^^y с номерами из множества
S^'9; матрица Fl^ имеет размеры v^x|Sl,e'9|, где v — размерность массива a, |Sl,e,9| —
число элементов множества s г ;
Li — число элементов массива а, помещающихся в строку кэша; предполагается, что строка массива помещается в кэш.
Лемма 1. Пусть для вхождения (l,p,q) массива а/ в правую часть оператора для всех S^9, 5<51,в,д, выполнены условия ^=1 (условие накладывается при наличии указанных 5) и v^|. Тогда в каждом потоке вычислений на разных итерациях циклов с параметрами j^, l^^1''3,9, используются разные элементы массива а.
Доказательство. Предположим противное утверждению леммы: на двух итерациях потока вычислений с различными параметрами циклов 5<5"в,д, используется один и тот же элемент данных. Это означает Fцi (¡J=Fцi на итерациях ,7 и ,7 потока вычислений с различными параметрами циклов £1,в,д; учтено, что т^=1 для
Так как матрица Fl ^ имеет |;~"в,д| линейно независимых столбцов, то можно выделить столько же линейно независимых строк. Матрицу, составленную из выделенных строк, обозначим F~. Обозначим через 7- и векторы размерности |£1'в,д|, построенные по векторам 7 и 7 выбором компонент с номерами из множества Тогда векторы 7-и различны и ,=Fl^ д 7^-Fl^ д 7^=0. Так как строки матрицы Fs составле-
ны из строк матрицы Fl^ , то выполняется равенство F~7~=F~J~. Из построения матрицы F~ следует, что она является невырожденной. Умножая последнее равенство слева на 1, получим 7~=7~ (противоречие). □
Отметим, что если F^p д^=0, "в,д, и итерации цикла с параметром выполняются параллельными потоками, то имеет место бродкаст.
Теорема 2. Пусть выполняются предположения леммы и не более чем один столбец F¿pqz имеет вид (0,^,0,р д^)Т, 0<\Ъ"^ . Если т^=1, цикл с параметром ^ вы-
полняется параллельно, циклы с параметрами в,д, F"^д ^Ф0, выполняются после-
довательно (с синхронизацией потоков вычислений между итерациями в случае т^>1), то независимо от выбора т^ 2 , в,д, совокупное число кэш-промахов по всем потокам вычислений является наименьшим. Сформулированные условия накладываются при наличии указанных столбцов.
Доказательство. Сделаем естественное предположение, что строка массива а^ выровнена в памяти по размеру строки кэша; в этом случае обращения к одной строке массива не уменьшат число кэш-промахов при обращении к данным, расположенным в начале следующей строки.
Условие леммы т^=1 для всех г^9, 5<51'в'g' в случае т^>1, цикл с параметром
^РлД
х^ =4Р,д ........ „„......__
цикла с параметром Так как выполняются предположения леммы, то в каждом по-
является последовательным и столбец F^^ линейно зависим со столбцами Fie г, Ses"3'9, можно заменить условием синхронизации потоков между итерациями
токе вычислений разные элементы данных используются на разных итерациях циклов с параметрами В каждом потоке вычислений на итерациях циклов с парамет-
рами (если 51'в,дфпр), используется один и тот же элемент данных, так как
столбцы F;вgr' являются нулевыми; итерации указанных циклов в потоке вы-
числений выполняются одна за другой. Следовательно, в каждом потоке на всех, кроме первой, итерациях, на которых используется элемент данных, он находится в кэше, то есть кэш-промахи, связанные с его использованием, отсутствуют. Если F"^q^=0, 5<51,в,д, и итерации цикла с параметром выполняются не одним, а параллельными потоками,
то осуществляется бродкаст и, следовательно, совокупное по потокам вычислений число кэш-промахов остается неизменным в полуварпе. Совокупное по потокам число кэш-
промахов также не изменяется при произвольном выборе г^2, 5е 5ФС, так как разные итерации параллельно выполняемых циклов , имеют доступ к разным строкам массива. Кроме того, циклы по ,5, выполняются последовательно (по
условию теоремы).
Рассмотрим цикл с параметром , (если имеется столбец указанного в формулировке теоремы вида). Использование элемента данных характеризуется пространственной локальностью относительно цикла с параметром так как используемые данные располагаются в одной строке массива на расстоянии Ьр^! элементов. Число используе-
мых строк кэша потоками равно
Ьг ,р ,9, и 20;,2
Ц
. Тогда в потоках вычислений имеется
Qi
ц
использований строк кэша без кэш-промахов при чтении элементов
строки массива (с учетом того, что строка массива а^ помещается в кэш и выполнено условие !Ь,в Их число наибольшее при наибольшем Q^2' то есть при г^=1.
Тогда совокупное по всем потокам вычислений число кэш-промахов, связанных с использованием элементов строки массива а, наименьшее (потребуется один обязательный кэш-промах для строк кэша). □
Теорема дает условия, при которых в потоках вычислений повторное использование элемента данных возникает только на последовательных итерациях алгоритма (следует из доказательства). Указывается выбор г^2, при котором число кэш-промахов наименьшее, а размеры г^2, для которых 5ег^9 или 5>5¿'в'^ могут быть заданы произвольно,
так как число кэш-промахов на этом вхождении будет одинаково по всем потокам вычислений. Выполнение условий теоремы свидетельствует о наличии пространственной локальности на итерациях цикла с параметром , (в различных потоках вычислений) и о
наличии временной локальности для осуществления бродкаста данного по параллельно выполняемым циклам с параметрами для которых Fцi д=0.
Пример 1 (продолжение). Рассмотрим блоки вычислений (рис. 3). Известно, что самый внешний цикл (цикл с параметром 1 Т^к) является последовательным, а внутренние циклы по ,2= и по 7з=7 могут выполняться параллельно. Имеем: F1 1 = 0 0 О ,
^13=1 0], К14=(1 0 0], 51Д,2=3, 51Д,3=2, 51,1,4=3. На вхождении (1,1,2) столбец , ,3 [ 10 0] , , [0 0 1 ,
Fl 1 2 1 является нулевым, поэтому, следуя теореме (предположения леммы справедливы и для теоремы) применительно к вхождению (1,в,с) =(1,1,2), положим ^ 2=1 (или положим Г1 2 максимально возможным — размеру блока ^ — и установим синхронизацию потоков между итерациями цикла по к). Далее на вхождении (1,1,2) имеем: Fl 1 22=(1 0)Т, 5=2, 5е г^9, поэтому ^ можно выбрать произвольно (но не больше максимально возможного размера блока); Fl 1 2з=(0 1)Т, поэтому выберем Гз2=1. На вхождении (1,1,3) имеем Fl 13 1=(0 1)Т, цикл по ^ выполняется последовательно, поэтому выбор Г1 2 не определяется теоремой; Fl 1 3 2=(1 0)Т, поэтому Г22 можно выбрать произволь-
но; так как ^1133=0, 5=3>51,1,3, то Т3 2 можно выбрать произвольно. На вхождении (1,1,4) имеем ^114 1=(1 0)Т, поэтому т12 можно выбрать произвольно; ^114 2=0, 5=2<51,1,4, поэтому выберем г22=1; F1 143=(0 1)Т, поэтому выберем г32=1.
Таким образом, суммируя сведения по каждому вхождению, выберем размеры потоков вычислений следующим образом: г^=1, ^=1, ^=1. Выбор размера потоков указанным образом определяет наименьшее в совокупности по потокам вычислений число кэш-промахов для каждого из вхождений (1,1,2), (1,1,3), (1,1,4) в отдельности. Бродкаст данных осуществляется на вхождении (1,1,3) для параллельного цикла по ^
(так как F11зз=0) и на вхождении (1,1,4) для параллельного цикла по ^2 (так как F1 1 4 2=0). Использование данных на итерациях цикла по ^ (в различных параллельных потоках вычислений) на вхождениях (1,1,2), (1,1,4) (F1 ! 2 з=F1 ! 4 3=(0 1)Т) характеризуется пространственной локальностью. Как уже отмечалось, можно положить т^ 2=^ и
установить синхронизацию потоков после каждой итерации к. Если потоки одного блока не имеют общих данных, то синхронизация не требуется.
4. Вычислительные эксперименты
Пример 1 (продолжение). В работе [10] реализован алгоритм Флойда—Уоршелла с 3D блоками. Как показано в подразделе 2.2, блочный алгоритм с 3D блоками обладает улучшенной, по сравнению с 2D-блочной версией, локальностью; это позволило уменьшить время реализации в 5-6 раз. В работе [11] алгоритм Флойда—Уоршелла с 3D блоками модифицирован, в том числе и с целью использования пространственной локальности данных; время реализации уменьшилось еще в 5 раз.
Пример 2. Рассмотрим алгоритм прямого хода метода Гаусса решения систем линейных алгебраических уравнений с использованием расширенной матрицы (рис. 4).
о к = 1, n - 1
do i = к + 1, n
do j = к + 1, n + 1
a(i,j) = a(ij) - (a(i,k) / a (k,k) ) * a(kj)
enddo
enddo
enddo
Рис. 4. Основная часть алгоритма Гаусса (без выбора ведущего элемента)
Исследование локальности этого алгоритма практически совпадает с приведенным в разделе 2 исследованием локальности алгоритма Флойда—Уоршелла. Некоторое несущественное дополнение к приведенным выкладкам связано с наличием ведущего элемента a(k,k).
В экспериментах использовалась видеокарта GeForce GT 860M со следующими характеристиками: объем глобальной памяти — 4295 МБ; максимальный размер разделяемой памяти в одном блоке — 49 КБ; количество 32-разрядных регистров в одном блоке — 65 536; количество потоков в варпе — 32; мультипроцессоров — 5; всего ядер — 640.
На рис. 5 продемонстрирована зависимость времени реализации алгоритма от ^ — параметра размера блока, характеризующего число записей в глобальную память. Чис-
п -1
ло записей результатов вычисления одного блока пропорционально Q1 =
поэтому
чем большие Гр тем меньше затраты времени на запись в глобальную память. Блоки имеют размер ^хт2хт3=т1 х64х64. В разделяемую память заносились только те элементы массива, которые использовались и для чтения, и для записи. Элементы, которые использовались только для чтения, читались из глобальной памяти. Пространственная локальность данных не использовалась.
Рис. 5. Зависимость времени реализации блочного алгоритма Гаусса от ^
Вычислительные эксперименты показали, что время реализации алгоритма уменьшается (до некоторого предела) с ростом Гр т.е. уменьшается при уменьшении числа
операций записи в глобальную память. Стабилизация времени связана с тем, что требуется определенное время на запуски ядра и на планирование потоков.
Заключение
Таким образом, в работе сформулированы и доказаны утверждения, позволяющие оценить объем коммуникационных операций, порождаемых разбиением множества итераций. Асимптотические оценки суммарного объема коммуникационных операций дают возможность ранжировать параметры размера блоков для минимизации объема коммуникационных операций. В работе также получены условия, при выполнении которых достигается наименьшее число кэш-промахов в потоках вычислений при наличии локальности данных.
Направления дальнейших исследований: разработка и программная реализация алгоритмов, позволяющих выбирать соотношения размеров блоков вычислений; получение условий и соотношений, точно характеризующих объем коммуникационных операций; обобщения достаточных условий наличия локальности данных, пригодные для применения в большем числе случаев; исследования, направленные на минимизацию объема коммуникационных операций между параллельными вычислительными процессами, реализуемыми на суперкомпьютерах с распределенной памятью; разработка и программная реализация параллельных алгоритмов для решения прикладных задач с использованием предлагаемого подхода.
Работа выполнена в рамках государственной программы научных исследований
Республики Беларусь «Конвергенция-2020», подпрограмма «Методы математического
моделирования сложных систем».
Литература
1. Лиходед Н.А. Полещук М.А. Метод ранжирования параметров размера блоков вычислений параллельного алгоритма / / Доклады НАН Беларуси. 2015. Т. 59, № 4. С. 25-33.
2. Kandemir M., Ramanujam J., Irwin M., Narayanan V., Kadayif I., Parikh A. A compiler based approach for dynamically managing scratch-pad memories in embedded systems / / IEEE Transactions on Computer-Aided Design. 2004. Vol. 23, No. 2. P. 243-260.
3. Baskaran M., Bondhugula U., Krishnamoorthy S., Ramanujam J., Rountev A., Sa-dayappan P. Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories / / Proceedings of the 13th ACM SIGPLAN Symposium on Principles and practice of parallel programming. Salt Lake City, USA, February 20-23, 2008.
4. Воеводин Вл.В., Воеводин Вад.В. Спасительная локальность суперкомпьютеров // Открытые системы. 2013. № 9. С. 12-15.
5. Воеводин В.В., Воеводин В.Вл. Параллельные вычисления. Санкт-Петербург: БХВ-Петербург, 2002. 608 с.
6. Xue J., Cai W. Time-minimal tiling when rise is larger than zero // Parallel Computing. 2002. Vol. 28, No. 5. P. 915-939.
7. Baskaran M., Ramanujam J., Sadayappan P. Automatic C-to-CUDA code generation for affine programs / / Proceedings of the Compiler Construction, 19th International Conference. Part of the Joint European Conferences on Theory and Practice of Software. Pa-phos, Cyprus, March 20-28, 2010.
8. Bondhugula U., Baskaran M., Krishnamoorthy S., Ramanujam J., Rountev A., Sadayappan P. Automatic transformations for communication-minimized parallelization and locality optimization in the polyhedral model // Lecture notes in computer science. 2008. No. 4959. P. 132-146.
9. Venkataraman G., Sahni S., Mukhopadhyaya S. A blocked all-pairs shortest-paths algorithm // J. Exp. Algorithmics. 2003. Vol. 8. P. 2.2.
10. Katz G.J., Kider J. All-pairs shortest-paths for large graphs on the GPU // Proceedings of the 23rd ACM SIGGRAPH/EUROGRAPHICS symposium on Graphics hardware. Sarajevo, Bosnia and Herzegovina: Eurographics Association. 2008. P. 47-55.
11. Lund B.D., Smith J.W. A multi-stage cuda kernel for floyd-warshall // CoRR abs/1001.4108. 2010.
DOI: 10.14529/cmse160307
ESTIMATE OF LOCALITY OF PARALLEL ALGORITHMS
IMPLEMENTED ON GPUS
© 2016 N.A. Likhoded, M.A. Paliashchuk
Belarusian State University (Nezavisimosti avenue 4, Minsk, 220030 Republic of Belarus)
E-mail: likhoded@bsu.by, poleschuma@bsu.by Received: 02.03.2016
The problem of obtaining blocks of operations and threads of parallel algorithm resulting in a smaller number of accesses to global memory and resulting in the efficient use of caches and shared memory graphics processor is investigated. We formulated and proved statements to assess the volume of communication transactions generated by alternative sizing of blocks, as well as to minimize the number of cache misses due to the use of temporal and spatial locality of data. The research is constructive and allows software implementation for practical use. Keywords: parallel computing, GPU, minimization of communications, temporal locality, spatial locality.
FOR CITATION
Likhoded N.A., Paliashchuk M.A. Estimate of Locality of Parallel Algorithms Implemented on GPUs. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2016. vol. 5, no. 3. pp. 96-111. (in Russian) DOI: 10.14529/cmse160307.
References
1. Likhoded N.A., Poleshchuk M.A. Method of Ranking Tiles Size Parameters of Parallel Algorithm. Doklady NAN Belarusi [Proceedings of the National Academy of Sciences of Belarus]. 2015. vol. 59, no. 4. pp. 25-33. (in Russian)
2. Kandemir M., Ramanujam J., Irwin M., Narayanan V., Kadayif I., Parikh A. A Compiler Based Approach for Dynamically Managing Scratch-Pad Memories in Embedded Systems. IEEE Transactions on Computer-Aided Design. 2004. vol. 23, no. 2. pp. 243-260. DOI: 10.1109/tcad.2003.822123.
3. Baskaran M., Bondhugula U., Krishnamoorthy S., Ramanujam J., Rountev A., Sa-dayappan P. Automatic Data Movement and Computation Mapping for Multi-Level Parallel Architectures with Explicitly Managed Memories. Proceedings of the 13th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. Salt Lake City, USA, February 20-23, 2008. DOI: 10.1145/1345206.1345210.
4. Voevodin Vl.V., Voevodin Vad.V. The Fortunate Locality of Supercomputers. Otkrytye sistemy [Open Systems]. 2013. no. 9. pp. 12-15. (in Russian)
5. Voevodin V.V., Voevodin Vl.V. Parallel'nye vychisleniya [Parallel Computing]. Sankt-Peterburg: BKhV-Peterburg, 2002. 608 p.
6. Xue J., Cai W. Time-Minimal Tiling when Rise is Larger than Zero. Parallel Computing. 2002. vol. 28, no. 5. pp. 915-939. DOI: 10.1016/s0167-8191(02)00098-4.
7. Baskaran M., Ramanujam J., Sadayappan P. Automatic C-to-CUDA Code Generation for Affine Programs. Proceedings of the Compiler Construction, 19th International Conference. Part of the Joint European Conferences on Theory and Practice of Software. Paphos, Cyprus, March 20-28, 2010. DOI: 10.1007/978-3-642-11970-5_14.
8. Bondhugula U., Baskaran M., Krishnamoorthy S., Ramanujam J., Rountev A., Sadayappan P. Automatic Transformations for Communication-Minimized Parallelization
H.A. JTnxo;i,<yi,. M.A. TTojieiuyK
and Locality Optimization in the Polyhedral Model. Lecture Notes in Computer Science. 2008. no. 4959. pp. 132-146. DOI: 10.1007/978-3-540-78791-4_9.
9. Venkataraman G., Sahni S., Mukhopadhyaya S. A Blocked All-Pairs Shortest-Paths Algorithm. J. Exp. Algorithmes. 2003. vol. 8. pp. 2.2. DOI: 10.1145/996546.996553.
10. Katz G.J., Kider J. All-Pairs Shortest-Paths for Large Graphs on the GPU. Proceedings of the 23rd ACM SIGGRAPH/EUROGRAPHICS Symposium on Graphics Hardware. Sarajevo, Bosnia and Herzegovina: Eurographics Association. 2008. pp. 47-55.
11. Lund B.D., Smith J.W. A Multi-Stage CUDA Kernel for Floyd-Warshall. CoRR abs/10014108. 2010.
