Исследуется задача получения блоков операций и потоков операций параллельного алгоритма, приводящих к меньшему числу обращений к глобальной памяти и к эффективному использованию параллельными потоками вычислений кэшей и разделяемой памяти графического процессора. Сформулированы и доказаны утверждения, позволяющие оценить объем коммуникационных операций, порождаемых альтернативными вариантами задания размеров блоков вычислений, а также минимизировать число промахов кэша за счет использования временной и пространственной локальности данных с учетом размера и длины строк кэша. Исследования конструктивны и допускают программную реализацию для практического использования. Ключевые слова: параллельные вычисления, графический процессор, минимизация объема коммуникационных операций, временная локальность, пространственная локальность.
Оценка локальности параллельных алгоритмов, реализуемых на графических процессорах* 2. Предварительные сведения
Будем считать, что алгоритм задан последовательной программой линейного класса1 [5]. Основную вычислительную часть такой программы составляют циклические конструкции; границы изменения параметров циклов задаются неоднородными формами, линейными по совокупности параметров циклов и внешних переменных. Предполагается, что в гнезде циклов имеется K выполняемых операторов Sβ и используется L массивов al размерностей νl. Область изменения параметров циклов (область итераций) для оператора Sβ и размерность этой области обозначим соответственно Vβ и nβ.
Вхождением (l,β,q) будем называть q-е вхождение массива al в оператор Sβ. Индексы элементов l-го массива, связанных с вхождением (l,β,q), выражаются функцией F lq вида
F lq (J ) FlqJ GlqN f lq J ( j1 jn ) V N Z s Flq Z ln Glq Z ls f lq Z l ,
где N Z s – вектор внешних переменных алгоритма, s – число внешних переменных.
Выполнение оператора Sβ при конкретных значениях β и вектора параметров цикла J будем называть операцией и обозначать Sβ(J). Пара вхождений (l,α,1) и (l,β,q) порождает истинную зависимость Sα(I)→Sβ(J), если: Sα(I) выполняется раньше Sβ(J); Sα(I) переопределяет (изменяет) элемент массива al, а Sβ(J) использует на вхождении (l,β,q) в качестве аргумента тот же элемент массива; между операциями Sα(I) и Sβ(J) этот элемент не переопределяется.
Зависимости (информационные связи) Sα(I)→Sβ(J) между операциями будем характеризовать векторами dα,β=J–I. Если размерности векторов J и I не совпадают, то J–I определяется как разность векторов меньшей из размерностей. Векторы dα,β часто называют векторами зависимостей: dα,β позволяет для операции Sβ(J) найти операцию Sα(I), от которой Sβ(J) зависит.
Пусть в гнезде циклов имеется наборов выполняемых операторов. Под набором операторов будем понимать один или несколько операторов, окруженных одним и тем же множеством циклов. Операторы и наборы операторов линейно упорядочены расположением их в записи алгоритма. Обозначим: V , 1 , – области изменения параметров циклов, окружающих наборы операторов, n – размерность области V , число циклов, окружающих -й набор операторов. Будем считать, что если оператор Sβ принадлежит набору операторов с номером β , то n =nβ.
Выделим блоки вычислений путем разбиения циклов, окружающих операторы. Пусть m J j1, jm2,,ijnn V j , M J j1, mj2,,ajnx V j , 1 n , – предельные значения изменения параметров циклов. Зададим размеры блоков вычислений натуральными числами r1 r . Параметр r n обозначает число значений параметра jζ, приходящихся на один блок -го набора операторов. Число частей Q , на которые при формировании блоков разбивается область значений параметра jζ цикла, окружающего -й набор операторов, находится согласно Q (M m 1) / r ( обозначает ближайшее «сверху» целое число). Нумеровать блоки вычислений будем по каждой координате в пределах от 0 до Q –1, 1 n .
Формализованным способом разбить множество операций алгоритма на блоки и потоки вычислений можно путем применения тайлинга – преобразования алгоритма для получения макроопераций [6, 7]. 1 В литературе используются также термины «аффинные гнезда циклов», «алгоритмы с аффинными зависимостями». Пример 1. Рассмотрим алгоритм Флойда-Уоршелла поиска кратчайших путей между всеми парами вершин графа (рис. 1).
do k = 1, n do i = 1, n do j = 1, n
a(i,j) = min(a(i,j), a(i,k) + a(k,j)) enddo enddo enddo
Рис. 2. Основная часть алгоритма Флойда-Уоршелла В гнезде циклов имеется один выполняемый оператор S1 (один набор операторов) и используется один массив a размерности 2. Область изменения параметров циклов (область итераций) V1=V 1 ={(k,i,j)Z3| 1≤k≤n, 1≤i≤n, 1≤j≤n} для оператора S1 имеет размерность 3. Для матриц Fl,β,q на вхождениях (l,β,q) имеем:
F1,1,1=F1,1,2= 0 1 0 , F1,1,3= 0 1 0 , F1,1,4= 1 0 0 .
0 0 1 1 0 0 0 0 1
В рассматриваемом примере векторы зависимостей d , будем для наглядности помечать
не номерами операторов, а элементами массивов, фигурирующими на порождающих
зависимости вхождениях. Например, вектор d a(i, j),a(i,k) порождается зависимостью между данными a(i,j)
на вхождении (l,α,1)=(
Выделим блоки вычислений. Напомним, блоки вычислений должны выполняться атомарно, независимо друг от друга. С учетом зависимостей алгоритма, такие блоки проще всего задать, если положить r1=1, r2 и r3 – параметры. Всего будет Q1×Q2×Q3 двумерных (2D) блоков, где Q1=n, Q2= n , Q3= n . Блочный алгоритм имеет следующий вид:
r2 r3 do k = 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 блоками 3. Оценка объема коммуникационных операций 3.1 Выражения для оценки объема операций чтения и записи Обозначим M max M m 1 – наибольшее число итераций циклов, участвующих в получении блоков. Для простоты записи будем использовать обозначение M без индекса и подразумевать набор операторов β при упоминании оператора Sβ. Отметим, что величина Qr имеет порядок M, поэтому r , Q могут быть величинами порядка M.
Определим термин «фиксированное данное массива» как конкретное, неизмененное содержимое соответствующей ячейки памяти. Введем в рассмотрение величины lq , ldq , коd торые определим следующим образом: пусть число фиксированных данных, используемых на d вхождении (l,β,q), оценивается величиной O(M lq ) , а число фиксированных данных, используемых на вхождении (l,β,q) при фиксированном значении цикла с параметром j оценивается величиной O(M ldq1) . Наличие индекса d в обозначениях подчеркивает, что если вхождение (l,β,q) порождает истинную зависимость Sα(I)→Sβ(J), то lq и ldq зависят от вектора dα,β. d Найти lq и ldq не представляет труда, если детально понимать реализуемый алгоритм.
d Формализованный способ получения величин lq , ldq изложен в работе [1]. Этот способ d использует функции F lq , а также функции, описывающие информационную структуру алгоритма – покрывающие функции графа алгоритма [5] (называемые еще h-преобразованиями [8]).
Отметим, что число фиксированных данных, используемых (в пространстве размерности nβ–1) на вхождении при фиксированном значении цикла с параметром jζ, по порядку либо равно ( ldq –1= ldq ), либо на 1 меньше ( ldq –1= ldq –1) числа фиксированных данных, используемых (в пространстве размерности nβ) при всех значениях цикла с параметром jζ.
Обозначим через d,,max максимальные по модулю значения компонент векторов dα,β. Сделаем естественное для практики предположение: компоненты векторов dα,β по модулю или не превосходят нескольких единиц, или неограниченно возрастают с ростом размера задачи. В первом случае d,,max также не превосходят нескольких единиц, во втором случае являются большими, значительно превосходящими r .
Теорема 1. Пусть вхождение (l,β,q) порождает истинную зависимость: определение некоторого данного происходит на вхождении (l,α,1) в левой части оператора Sα, а использование – на вхождении (l,β,q) в правой части оператора Sβ, причем в окружении операторов Sα и Sβ имеется цикл с параметром jζ.
Тогда, при получении блоков вычислений для реализации алгоритма на графическом процессоре, разбиение итераций цикла jζ порождает коммуникационные операции чтения и записи, объем которых имеет следующие оценки:
d O(QM lq1) операций чтения и операций записи, если 0 d,,max r ;
d d O(QM lq ) операций чтения и O(M lq ) операций записи, если d ,,max r , ldq ≠ ldq ;
d O(M lq ) операций чтения и операций записи, если d,,max r , ldq = ldq ; не требуется ни операций чтения, ни операций записи, если d,,max =0. В случае, когда вхождение (l,β,q) не порождает истинной зависимости (происходит обращение к входным данным) или цикл с параметром jζ имеется в окружении только оператора Sβ, оценки следующие:
d O(QM lq ) операций чтения, если ldq ≠ ldq ;
d
O(M lq ) операций чтения, если ldq = ldq .
Доказательство. Утверждения теоремы оценивают объем коммуникационных операций, порождаемых разбиением итераций цикла jζ, в зависимости от значений ζ-й координаты вектора dα,β и величин ldq , lq . Напомним, множество итераций цикла с параметром jζ разбиваетd ся на Q частей, каждая часть имеет некоторый номер jgl ( 0 jgl Q 1 ) и содержит r итераций (кроме, возможно, части с номером Q 1).
Пусть 0 d,,max r , определение данного происходит при выполнении операции S(I (i1in )) , а использование – при выполнении операции S(J ( j1 jn )) .
Число фиксированных данных, используемых на вхождении (l,β,q) при фиксированном значении цикла с параметром jζ (т.е. в пространстве размерности nβ–1), оценивается величиной O(M ldq1) . Тогда число фиксированных данных, используемых на вхождении (l,β,q) при всех значениях цикла с параметром jζ (в пространстве размерности nβ), оценивается величиной
M d ,,max O(M ldq1) O(M ldq ) (напомним, величина d ,,max не превышает нескольких единиц). d С другой стороны, оценка числа таких данных есть O(M lq ) . Таким образом, в рассматриваемом случае верно ldq = ldq .
В каждой из Q частей итераций цикла jζ число фиксированных данных, которые определяются при одном значении jζgl , но используются при другом, оцениваются величиной d ,,maxO(M ldq 1) , причем независимо от значения r величина d ,,max , как было условлено, не превышает нескольких единиц. Суммарный объем коммуникационных операций чтения и операций записи по всем Q частям определяется оценкой QO(M ldq1) = O(QM lq1) . d Пусть d ,,max r . Если ldq ≠ lq (т.е. ldq –1= ldq ), то число используемых на вхожd дении (l,β,q) фиксированных данных в каждой части и во всех Q частях оценивается одинакоd вой величиной O(M lq ) . Этой величиной следует оценить объем коммуникационных операций записи. Оценка объема коммуникационных операций чтения по всем Q частям есть d d d QO(M lq ) = O(QM lq ) . Если ldq = lq , то число используемых фиксированных данных в каждой из Q частей оценивается величиной O(rM ldq1) . Суммарный объем коммуникационных операций чтения и операций записи по всем Q частям определяется оценкой QO(rM ldq1) = O(QrM lq1) = O(M lq ) .
d d Если d,,max =0, то любое фиксированное данное определяется и используется при одном и том же значении параметра цикла jζgl . Коммуникационных операций не требуется.
Пусть вхождение (l,β,q) не порождает истинной зависимости (происходит обращение к входным данным) или порождает истинную зависимость, но цикл с параметром jζ имеется в окружении только оператора Sβ. Если ldq ≠ lq (т.е. ldq –1= ldq ), то объем коммуникациd на вхождении при фиксированном значении одного из циклов с параметрами j1=i, j2=j, j3=k, оценивается величинами O(M 1d121) =O(n2).
Для третьего вхождения a(i,k) (использование при фиксированном k на итерациях (i,j,k) элементов одного столбца) d a(i, j),a(i,k) =(1,0,j–k), 113 =2, 1d1,13 =2, 1d1,23 =2, 1d1,33 =3. d Для четвертого вхождения a(k,j) (использование при фиксированном k на итерациях (i,j,k) элементов одной строки) d a(i, j),a(k, j) =(1,i–k,0), 114 =2, 1d1,14 =2, 1d1,24 =3, 1d1,34 =2. d Оценим, используя теорему 1, объем коммуникационных операций, порождаемых разбиением итераций циклов.
Пусть ζ=1. Тогда для всех вхождений d,,max = d, =1. Имеет место случай
d
0 d,,max r , требуется O(QM lq1) операций чтения и операций записи: для второго
вхождения O(Q1n2) операций, для третьего и четвертого вхождений – O(Q1n) операций.
Пусть ζ=2. Для второго и третьего вхождений d ,,max = d, =0, поэтому не требуется ни
операций чтения, ни операций записи. Для четвертого вхождения (случай d ,,max r , ldq ≠
d d d
lq , (l,β,q)=(
Пусть ζ=3. Для второго и четвертого вхождений d ,,max = d, =0, не требуется ни
опера
ций чтения, ни операций записи. Для третьего вхождения (случай d ,,max r , ldq ≠ ldq ,
(l,β,q)=(
Утверждения теоремы 1 позволяют найти асимптотику суммарного объема коммуникационных операций, порождаемых разбиением множества итераций j на Q частей, и, следовательно, дают возможность ранжировать (устанавливать соотношение размеров относительно друг друга) параметры размера блоков вычислений для минимизации объема коммуникационных операций. При ранжировании параметров размера блоков вычислений параллельного алгоритма, реализуемого на графическом процессоре, следует в оценках объема коммуникационных операций принимать во внимание только оценки, содержащие множитель Q ; следует также учитывать, что величины r и Q связаны обратно пропорциональной зависимостью. Пример 1 (продолжение). Для ранжирования параметров r1, r2 и r3 размера блоков вычислений, выполняемых на графическом процессоре, укажем, следуя результатам предыдущего раздела, оценки объема коммуникационных операций, содержащие множитель Q – число частей, на которые при формировании блоков разбивается область значений параметра jζ.
Если ζ=1, то требуется O(Q1n2) операций чтения и операций записи для второго вхождения и O(Q1n) операций чтения и операций записи для третьего и четвертого вхождений. Если ζ=2, то требуется O(Q2n2) операций чтения для четвертого вхождения. Если ζ=3, то требуется O(Q3n2) операций чтения для третьего вхождения.
Таким образом, сравнивая оценки коммуникационных издержек и по чтению и по записи, приходим к выводу, что параметр размера блока вычислений r1 увеличивать более выгодно, чем параметры r2 и r3. В то же время, наиболее простой блочный алгоритм Флойда-Уоршелла (рис. 2) содержит 2D блоки, для которых r1=1. Блочный алгоритм с 3D блоками (рис. 3), для которых r1=r2=r3=r>1 получен в работе [9]. Представляет также интерес разработка блочного алгоритма с 3D блоками, для которых r1>r2, r1>r3.
do k = 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 блоки блочного алгоритма Флойда-Уоршелла 4. Условия наличия локальности данных
Потоки одного блока выполняются на мультипроцессоре частями-пулами, называемыми варп. Варп содержит до 32 потоков (два полуварпа по 16 потоков). Если несколько потоков одного полуварпа обращаются к одному и тому же банку разделяемой памяти, константной памяти или кэша текстурной памяти для доступа к различным данным, то происходит конфликт. Но конфликта не происходит, если несколько потоков одного полуварпа обращаются к одному и тому же данному (broadcast). В этом случае запрашиваемое данное передается только один раз, поэтому трафик сокращается в 16 раз (для архитектуры Fermi в 32 раза, так как объединение запросов в память происходит для 32 потоков).
Текстурная память может быть использована эффективно, если каждый поток полуварпа запрашивает близко расположенные в памяти данные, то есть если доступ к памяти характеризуется пространственной локальностью. Наличие пространственной локальности позволяет эффективно использовать кэш текстурной памяти и разделяемую память параллельными потоками.
Зададим в блоках вычислений потоки вычислений посредством выделения блоков второго уровня. Размеры блоков второго уровня зададим натуральными числами r1,2, r2,2, …, rnβ,2. Параметр rξ,2 обозначает число значений параметра jξ, приходящихся на один поток (на один блок второго уровня). Обозначим через Qξ,2 число частей, на которые при формировании потоков разбивается область значений цикла с параметром jξ, приходящихся на один блок; Qξ,2= r / r,2 . Выбор размеров rξ,2, оптимальных по числу кэш-промахов, определяется в значительной степени наличием временной локальности и пространственной локальности данных. Введем обозначения: Fl,β,q,ξ – столбец с номером ξ матрицы Fl,β,q; ξl,β,q – номер ненулевого столбца матрицы Fl,β,q, правее которого находятся только нулевые столбцы матрицы; если самый правый столбец Fl,β,q,nβ не нулевой, то по определению ξl,β,q=nβ; Ξl,β,q – множество номеров линейно независимых столбцов матрицы Fl,β,q (если таких множеств более одного, то зафиксировать любое);
Flq – матрица, составленная из столбцов матрицы Fl,β,q с номерами из множества Ξl,β,q; матрица Flq имеет размеры νl×|Ξl,β,q|, где νl – размерность массива al, |Ξl,β,q| – число элементов множества Ξl,β,q;
Ll – число элементов массива al, помещающихся в строку кэша; предполагается, что строка массива помещается в кэш.
Лемма. Пусть для вхождения (l,β,q) массива al в правую часть оператора Sβ для всех ξ Ξl,β,q, ξ≤ξl,β,q, выполнены условия rξ,2=1 (условие накладывается при наличии указанных ξ) и νl≥|Ξl,β,q|. Тогда в каждом потоке вычислений на разных итерациях циклов с параметрами jξ, 1≤ξ≤ξl,β,q, используются разные элементы массива al.
Доказательство. Предположим противное утверждению леммы: на двух итерациях потока вычислений с различными параметрами циклов jξ, ξ≤ξl,β,q, используется один и тот же элемент данных. Это означает Fl,β,qJ=Fl,β,qĴ на итерациях J и Ĵ потока вычислений с различными параметрами циклов jξ, ξΞl,β,q; учтено, что rξ,2=1 для ξΞl,β,q, ξ≤ξl,β,q.
Так как матрица Flq имеет |Ξl,β,q| линейно независимых столбцов, то можно выделить столько же линейно независимых строк. Матрицу, составленную из выделенных строк, обозначим FΞ. Обозначим через JΞ и ĴΞ векторы размерности |Ξl,β,q|, построенные по векторам J и Ĵ выбором компонент с номерами из множества Ξl,β,q. Тогда векторы JΞ и ĴΞ различны и Fl,β,qJ– Fl,β,qĴ= Flq JΞ– Flq ĴΞ=0. Так как строки матрицы FΞ составлены из строк матрицы Flq , то выполняется равенство FΞJΞ=FΞĴΞ. Из построения матрицы FΞ следует, что она является невырожденной. Умножая последнее равенство слева на (FΞ)–1, получим JΞ=ĴΞ (противоречие). □ Отметим, что если Fl,β,q,ξ=0, ξ<ξl,β,q, и итерации цикла с параметром jξ выполняются параллельными потоками, то имеет место бродкаст.
Теорема 2. Пусть выполняются предположения леммы и не более чем один столбец Fl,β,q,ζ имеет вид (0,…,0,bl,β,q,ζ)T, 0<|bl,β,q,ζ|rζ,2<Ll. Если rζ,2=1, цикл с параметром jζ выполняется параллельно, циклы с параметрами jξ, ξΞl,β,q, Fl,β,q,ξ≠0, выполняются последовательно (с синхронизацией потоков вычислений между итерациями в случае rξ>1), то независимо от выбора rξ,2, ξ>ξl,β,q, совокупное число кэш-промахов по всем потокам вычислений является наименьшим. Сформулированные условия накладываются при наличии указанных столбцов.
Доказательство. Сделаем естественное предположение, что строка массива al выровнена в памяти по размеру строки кэша; в этом случае обращения к одной строке массива не уменьшат число кэш-промахов при обращении к данным, расположенным в начале следующей строки.
Условие леммы rξ,2=1 для всех ξΞl,β,q, ξ≤ξl,β,q, в случае rξ>1, цикл с параметром jξ, ξ≤ξl,β,q, является последовательным и столбец Fl,β,q,ξ линейно зависим со столбцами Fl,β,q,ξ, ξ Ξl,β,q, можно заменить условием синхронизации потоков между итерациями цикла с параметром jξ. Так как выполняются предположения леммы, то в каждом потоке вычислений разные элементы данных используются на разных итерациях циклов с параметрами jξ, ξ≤ξl,β,q. В каждом потоке вычислений на итерациях циклов с параметрами jξ, ξ>ξl,β,q (если ξl,β,q≠nβ), используется один и тот же элемент данных, так как столбцы Fl,β,q,ξ, ξ>ξl,β,q, являются нулевыми; итерации указанных циклов в потоке вычислений выполняются одна за другой. Следовательно, в каждом потоке на всех, кроме первой, итерациях, на которых используется элемент данных, он находится в кэше, то есть кэш-промахи, связанные с его использованием, отсутствуют. Если Fl,β,q,ξ=0, ξ<ξl,β,q, и итерации цикла с параметром jξ выполняются не одним, а параллельными потоками, то осуществляется бродкаст и, следовательно, совокупное по потокам вычислений число кэшпромахов остается неизменным в полуварпе. Совокупное по потокам число кэш-промахов также не изменяется при произвольном выборе rξ,2, ξ Ξl,β,q, ξ≠ζ, так как разные итерации параллельно выполняемых циклов jξ имеют доступ к разным строкам массива. Кроме того, циклы по jξ, ξΞl,β,q, Fl,β,q,ξ≠0, выполняются последовательно (по условию теоремы).
Рассмотрим цикл с параметром jζ (если имеется столбец указанного в формулировке теоремы вида). Использование элемента данных характеризуется пространственной локальностью относительно цикла с параметром jζ, так как используемые данные располагаются в одной строке массива на расстоянии |bl,β,q,ζ| элементов. Число используемых строк кэша потоками равно blβq,ζ rζ,2Qζ,2 . Тогда в потоках вычислений имеется Qζ,2– blβq,ζ rζ,2Qζ,2 использований Ll Ll строк кэша без кэш-промахов при чтении элементов строки массива (с учетом того, что строка массива al помещается в кэш и выполнено условие |bl,β,q,ζ|rζ,2<Ll). Их число наибольшее при наибольшем Qζ,2, то есть при rζ,2=1. Тогда совокупное по всем потокам вычислений число кэшпромахов, связанных с использованием элементов строки массива al, наименьшее (потребуется один обязательный кэш-промах для строк кэша). □ Теорема дает условия, при которых в потоках вычислений повторное использование элемента данных возникает только на последовательных итерациях алгоритма (следует из доказательства). Указывается выбор rξ,2, при котором число кэш-промахов наименьшее, а размеры rξ,2, для которых ξΞl,β,q или ξ>ξl,β,q, могут быть заданы произвольно, так как число кэшпромахов на этом вхождении будет одинаково по всем потокам вычислений. Выполнение условий теоремы свидетельствует о наличии пространственной локальности на итерациях цикла с параметром jζ (в различных потоках вычислений) и о наличии временной локальности для осуществления бродкаста данного по параллельно выполняемым циклам с параметрами jξ, для которых Fl,β,q,ξ=0.
Пример 1 (продолжение). Рассмотрим блоки вычислений (рис. 3). Известно, что самый
внешний цикл (цикл с параметром j1, j1=k) является последовательным, а внутренние циклы по
j2=i и по j3=j могут выполняться параллельно. Имеем: F1,1,2= 0 1 0 , F1,1,3= 0 1 0 , F1,1,4=
0 0 1 1 0 0
1 0 0 , ξ1,1,2=3, ξ1,1,3=2, ξ1,1,4=3. На вхождении (
Таким образом, суммируя сведения по каждому вхождению, выберем размеры потоков
вычислений следующим образом: r1,2=1, r2,2=1, r3,2=1. Выбор размера потоков указанным образом
определяет наименьшее в совокупности по потокам вычислений число кэш-промахов для
каждого из вхождений (
Пример 1 (продолжение). В работе [10] реализован алгоритм Флойда-Уоршелла с 3D блоками. Как показано в подразделе 3.2, блочный алгоритм с 3D блоками обладает улучшенной, по сравнению с 2D-блочной версией, локальностью; это позволило уменьшить время реализации в 5-6 раз. В работе [11] алгоритм Флойда-Уоршелла с 3D блоками модифицирован, в том числе и с целью использования пространственной локальности данных; время реализации уменьшилось еще в 5 раз.
Пример 2. Рассмотрим алгоритм прямого хода метода Гаусса решения систем линейных алгебраических уравнений с использованием расширенной матрицы (рис. 4).
do k = 1, n – 1 do i = k + 1, n do j = k + 1, n + 1
a(i,j) = a(i,j) – (a(i,k) / a(k,k)) * a(k,j) enddo enddo enddo
Рис. 4. Основная часть алгоритма Гаусса (без выбора ведущего элемента) Исследование локальности этого алгоритма практически совпадает с приведенным в разделе 3 исследованием локальности алгоритма Флойда-Уоршелла. Некоторое несущественное дополнение к приведенным выкладкам связано с наличием ведущего элемента a(k,k).
В экспериментах использовалась видеокарта GeForce GT 860M со следующими характеристиками: объем глобальной памяти – 4295 МБ; максимальный размер разделяемой памяти в одном блоке – 49 КБ; количество 32-разрядных регистров в одном блоке – 65 536; количество потоков в варпе – 32; мультипроцессоров – 5; всего ядер – 640.
На рис. 5 продемонстрирована зависимость времени реализации алгоритма от r1 – параметра размера блока, характеризующего число записей в глобальную память. Число записей ре n 1 зультатов вычисления одного блока пропорционально Q1 r1 . Использование 3D блоков размером r1×r2×r3=r1×64×64 с большим значением r1 позволяет уменьшить затраты времени на доступ в глобальную память. В разделяемую память заносились только те элементы массива, которые использовались и для чтения, и для записи. Элементы, которые использовались только для чтения, читались из глобальной памяти. Пространственная локальность данных не использовалась.
Рис. 5. Зависимость времени реализации блочного алгоритма Гаусса от r1 в глобальную память. Стабилизация времени связана с тем, что требуется определенное время на запуски ядра и на планирование потоков. 6. Заключение
Таким образом, в работе сформулированы и доказаны утверждения, позволяющие оценить объем коммуникационных операций, порождаемых разбиением множества итераций. Асимптотические оценки суммарного объема коммуникационных операций дают возможность ранжировать параметры размера блоков для минимизации объема коммуникационных операций. В работе также получены условия, при выполнении которых достигается наименьшее число кэшпромахов в потоках вычислений при наличии локальности данных.
Направления дальнейших исследований: разработка и программная реализация алгоритмов, позволяющих выбирать соотношения размеров блоков вычислений; получение условий и соотношений, точно характеризующих объем коммуникационных операций; обобщения условий наличия локальности данных, пригодные для применения в большем числе случаев; исследования, направленные на минимизацию объема коммуникационных операций между параллельными вычислительными процессами, реализуемыми на суперкомпьютерах с распределенной памятью; разработка и программная реализация параллельных алгоритмов для решения прикладных задач с использованием предлагаемого подхода. Литература 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. 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. 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. 11. Lund B.D., Smith J.W. A multi-stage cuda kernel for floyd-warshall // CoRR abs/1001.4108. 2010. Estimate of locality of parallel algorithms implemented on GPUs
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. 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.