CC BY
16
УДК 004.021, 519.685 DOI: 10.14529/cmse200105
АНАЛИТИЧЕСКОЕ МОДЕЛИРОВАНИЕ МАТРИЧНО-ВЕКТОРНОГО ПРОИЗВЕДЕНИЯ НА МНОГОЯДЕРНЫХ ПРОЦЕССОРАХ*
© 2020 Е.Н. Акимова1,2, Р.А. Гареев1
1 Уральский федеральный университет им. Б.Н. Ельцина (620002 Екатеринбург, у л. Мира, д. 19),
2Институт математики и механики им. Н.Н. Красовского УрО РАН (620990 Екатеринбург, ул. Софьи Ковалевской, д. 16)
E-mail: aenl5@yandex.ru, gareevroman@gmail.com Поступила в редакцию: 31.01.2020
Эффективная реализация матрично-векторного произведения имеет существенную практическую значимость в областях машинного обучения, интеллектуального анализа данных, квантовой химии, математической физики, численных методов линейной алгебры, высокопроизводительных вычислений и др. В данной работе представлен алгоритм автоматизированной оптимизации матрично-векторного произведения по времени выполнения, использующийся на этапе компиляции без ручной настройки и автонастройки. Алгоритм основан на моделировании вычислений на гипотетическом многоядерном процессоре, предложенном авторами, с применением полиэдрального представления. В отличие от подходов, основанных на ручной настройке и автонастройке, алгоритм может применяться для создания новых оптимизированных реализаций матрично-векторного произведения в условиях недоступности целевой архитектуры и ограниченности времени выполнения. Алгоритм использован для оптимизации программного кода, реализующего решение структурной обратной задачи гравиметрии о нахождении поверхности раздела сред методом Левенберга— Марквардта. Проведено сравнение производительности полученной реализации с реализациями на основе оптимизированных библиотек линейной алгебры Intel MKL, BLIS, OpenBLAS. Результаты численных экспериментов показывают сравнимость предложенного алгоритма по эффективности с подходами, созданными с использованием ручной настройки при доступе к целевым архитектурам процессоров.
Ключевые слова: компиляторы, линейная алгебра, матрично-векторные операции, аналитическое моделирование, обратная задача гравиметрии.
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Акимова Е.Н., Гареев Р.А. Аналитическое моделирование матрично-векторного произведения на многоядерных процессорах // Вестник ЮУрГУ. Серия: Вычислительная математика и информатика. 2020. Т. 9, № 1. С. 69-82. DOI: 10.14529/cmse200105.
Введение
Матрично-векторное произведение (matrix-vector multiplication, MVM) и матричное произведение (matrix-matrix multiplication, МММ) широко используются для решения прикладных задач. В частности, MVM вычисляется в процессе обрезки весов (weight pruning) в глубинных нейронных сетях (Deep Neural Networks, DNN), используемых в машинном обучении [1]. MVM применяется в интеллектуальном анализе графов [2]. MVM вычисляется при определении коэффициента скорости реакции в квантовой химии [3]. MVM и МММ используются для решения широкого ряда задач математической физики и, в частности, математической геофизики, численными методами [4]. Минимизация времени выполнения MVM является отдельным предметом изучения высокопроизводительных вычислений [5, 6].
‘Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии (ПаВТ) 2020»
Современные подходы, используемые для оптимизации MVM по времени выполнения, основаны на эмпирическом поиске. Их можно разделить на два вида: ручная настройка (manual-tuning) и автонастройка (auto-tuning) [7]. Для выполнения ручной настройки для отдельно рассматриваемой архитектуры компьютера специалист, обладающий знаниями особенностей используемых аппаратных средств, создает новую реализацию требуемой операции. В случае автонастройки для каждой отдельно рассматриваемой компьютерной архитектуры выполняется перебор значений параметров заранее созданной реализации для нахождения их оптимальных значений. Это делает невозможным их применение в процессе кросс-компиляции и в других условиях недоступности целевой архитектуры процессора. Выполнение ручной настройки и автонастройки может потребовать значительных временных затрат. Вследствие этого ручная настройка и автонастройка не могут использоваться, когда время выполнения ограничено. Таким примером являются оптимизации, применяемые компилятором по умолчанию. С целью возможности их частого использования на этапе разработки и компиляции программ предпочтительным являются подходы, требующие наименьшего времени выполнения.
Результаты [7] в области моделирования МММ на гипотетическом процессоре показывают, что эмпирический поиск необязателен для достижения высокой производительности. Рассматриваемые подходы для оптимизации МММ не могут быть использованы для MVM из-за меньшей доли вычислений и большей доли обращений к памяти [8].
В работах [9, 10] представлены алгоритмы оптимизации по времени выполнения обобщенных МММ и сверток тензоров (tensor contraction, ТС) в процессе компиляции.
В данной работе предложен алгоритм оптимизации по времени выполнения MVM без эмпирического поиска. Оптимизация выполняется с использованием полиэдрального представления. Для определения значений параметров оптимизации применяется моделирование вычисления MVM на гипотетическом многоядерном процессоре. Алгоритм использован для оптимизации программы, реализующей решение обратной структурной задачи гравиметрии о нахождении поверхности раздела сред методом Левенберга—Марквардта.
Полученный программный код сравнивается по производительности с программными кодами, оптимизированными с помощью библиотек линейной алгебры, созданных с использованием ручной настройки при доступе к целевым архитектурам процессоров (Intel MKL [11], BLIS [12], OpenBLAS [13]). Преимущество предложенного алгоритма состоит в возможности его применения для создания новых оптимизированных реализаций MVM в условиях недоступности целевой архитектуры и ограниченности времени выполнения.
Численные эксперименты проводились на узле с двумя 18-ядерными процессорами Intel Xeon Е5-2697 v4 суперкомпьютера «Уран» (Институт математики и механики им. Н.Н. Красовского УрО РАН, Екатеринбург, Россия). Показана сравнимость нашего подхода по производительности с подходами, основанными на эмпирическом поиске.
Статья организована следующим образом. Раздел 1 посвящен описанию алгоритма автоматизированной оптимизации матрично-векторного произведения. В разделе 2 выполнена постановка структурной обратной задачи гравиметрии, описан алгоритм её решения, приведены результаты сравнения производительности полученной реализации решения задачи гравиметрии с реализациями на основе оптимизированных библиотек линейной алгебры. В заключении содержатся выводы и направления дальнейших исследований.
1. Автоматизированная оптимизация MVM
В данном разделе описывается алгоритм автоматизированной оптимизации (ААО) MVM, использующий аналитическое моделирование и оптимизацию полиэдрального представления. Описываемые методы могут использоваться во время кросс-компиляции, выполняемой промышленными компиляторами, и в других случаях, когда время выполнения компиляции ограничено либо архитектура целевого процессора недоступна.
Ручная оптимизация широко применяется для MVM и других операций линейной алгебры, выполняемых экспертами в области ее приложений [5]. Для того чтобы упростить разработку новых высокопроизводительных реализаций под новые архитектуры процессоров, специалистами была предложена автонастройка и аналитическое моделирование [7]. Данные подходы используются для нахождения значений параметров созданных реализаций, обеспечивающих наилучшую производительность программ.
В случае автонастройки время выполнения рассматриваемой реализации измеряется для каждого набора значений параметров. Аналитическое моделирование позволяет смоделировать потребление ресурсов гипотетического процессора (ГП) [7] с целью установления зависимостей между производительностью и значениями параметров рассматриваемой реализации. Для определения наилучших значений параметров реализации MVM в работе предложена модификация аналитической модели, которая использовалась в библиотеке BLIS для МММ.
1.1. Модель гипотетического процессора
ГП описывается следующим образом [7].
• Архитектура загрузки/сохранения и векторные регистры: данные должны быть загружены в регистры процессора перед тем как с ними могут быть выполнены вычисления. Имеется Адес векторных регистров. Каждый из них может содержать Ауес значений размера SbАТА- Если Адес равно 0, Avec присваивается значение 1. Предполагается, что команды для работы с памятью могут выполняться одновременно с командами арифметики с плавающей точкой.
• Векторные инструкции: пропускная способность процессора составляет АудмА векторных операций смешанного сложения и умножения (vector fused multiply-add instruction, VFMA) за один такт. Отдельная VFMA выполняет Avec умножений и сложений, составляющих VFMA. Минимальное количество тактов ГудмА) которое должно быть совершено перед началом выполнения новой зависимой по данным VFMA-инструкции, определяет задержку VFMA-инструкции.
• Кэш-память: весь кэш данных — это множественно-ассоциативный кэш с политикой вытеснения последнего по времени использования (Least Recently Used, LRU). Каждый уровень кэша L* характеризуется следующими параметрами: размер линии кэша Сд, степень ассоциативности Шд, количество множеств Ад, размер 5д, где 5д = АдСдИ7Li- В случае полностью ассоциативного кэша Ад = 1.
Авторами предложена модификация ГП, которая определяет инструкции предвыборки, помещающие данные в кэш первого уровня.
• Инструкции предвыборки: за один такт процессора может быть выполнено
-^prefetch инструкций. Задержка каждой инструкции составляет -Lprefetch тактов процессора. Каждая инструкция может загрузить данные, размер которых равен Cbt-
Помимо инструкций предвыборки модификация определяет минимальное количество тактов -LvloaDj которое должно быть совершено перед началом выполнения новой зависимой по данным векторной инструкции загрузки данных на векторный регистр процессора.
Вычислительная сложность MVM составляет 0(N2), что в N раз меньше, чем вычислительная сложность МММ. Вследствие этого каждый элемент умножаемой матрицы используется только 0(1) раз, а не 0(N) как в случае МММ. Следовательно, стоимость загрузки элементов матрицы из памяти на регистры процессора преобладает над количеством выполняемых операций. Предвыборка данных, загружающая элементы матрицы в кэш первого уровня до момента их использования, может уменьшить количество промахов кэша и, следовательно, уменьшить стоимость загрузки элементов матрицы. Как следствие, эффективная предвыборка необходима в случае MVM.
1.2. Алгоритмы вычисления MVM
Алгоритм 1 вычисления MVM для случая транспонированной матрицы имеет следующий вид:
Алгоритм 1: Вычисление MVM. Случай транспонированной матрицы
i begin
2
3
4
5
6
7
8 9
10
11
12
13
14
15
16
for jc = 0...N step Nc do for ic = 0...M step Mc do for jb = 0...1VC step Щ do
Выполняем предвыборку Mc x Щ элементов матрицы А с шагом D for ib — О...Мс step 1 do I — ic + ib acc(0:Nr — 1) = x(I) for jr — О...Щ step Nr do
J = jc + jb + jr ■ jc + jb + jr + Aj. — 1
y(JT) += A(/,J) • acc(0 : Nr — 1) end
end
end
end
end
17 end
Алгоритм 2 вычисления MVM для случая нетранспонированной матрицы имеет следующий вид:
Алгоритм 2: Вычисление MVM. Случай нетранспонированной матрицы
i begin
2
3
4
5
6
7
8 9
10
11
12
13
14
15
16
17
18
for jc = 0...N step Nc do for ic = 0...M step Mc do for jr — O...Nc step Nr do
if jr mod oCli == 0 then
J о ПАТА
4 CLl
Предвыборка Mc x gDA^A элементов А с шагом D 4 Cl
Предвыборка gDA^A элементов x с шагом D end end
JS — jc jr • jc A jr “b Nr 1 acc(0,0 : Nr — 1) +— A(ic,S) ■ ж(Л)
acc(Mc — 1, 0 : Nr — 1) +— A(ic + Mc — 1, J) • ж(Л) end
y{ic) += EiIro_1 acc(°> *)
y(ic + MC- 1) += E£о 1 acc(Mc - 1, i)
end
19 end
В процессе выполнения алгоритмов 1 и 2 матрицы и векторы разбиваются на блоки. Это влияет на выполняемую предвыборку данных, загружающую элементы матрицы в кэш первого уровня, а также позволяет использовать элементы блоков, хранящихся в кэш-памяти, повторно. Размеры блоков и размер шага предвыборки являются параметрами алгоритма. Их оптимальные значения определяются с помощью гипотетического процессора, описанного в разделе 1.1.
В случае транспонированной матрицы формулы для нахождения значений параметров имеют следующий вид:
Щ = min
| max {AVfma-Lvfma-Nvec, {Nreg - 2)AfyEc},
AVec5data
Щес
Nr = NyEC, Mc = min { Wl2 - 1, WLl}, Nc = Nb2 Cl2/Sdata ,
jj _ -^prefetch
\ Nlc \ Nb^DATA/(-'Ll 1 /Nprefetch] + Mc(\ Nb/(NyEcNyFMA)] + ^VLOAd)
В случае нетранспонированной матрицы формулы для нахождения значений параметров имеют следующий вид:
Nc = NL2Cl2/Sdata, Nr = Nyeg \j/mm{WLl,WL2 - l}] , Mc = frNym/N^
7 = -^YFMA^YFMA,
-^prefetch
[4(Мс + l)/-/Vprefetch] + 4Сщ ( \McNr/(-/VvFMA-^VvEc)! + AvLOAd)/(-^ВАТА-^г)
Значение параметров -/Vreg, AfyEc, -Sdata, AfyFMA, Tvfma, АД., Cl<, 1ТД, ^prefetch, -^prefetch) -^VLOAD являются параметрами модели гипотетического процессора, описанной в разделе 1.1.
В реальных случаях параметры определяются техническими характеристиками целевого процессора.
1.3. Комплекс программ
Для автоматизированного распознавания и оптимизации MVM в промышленных компиляторах применяется полиэдральная компиляция. Основой такого подхода служит полиэдральное представление программы [14].
Определение. Полиэдральная модель — это математическое представление (фреймворк) для моделирования и оптимизации доступа к памяти, производимого в группах вложенных циклов, не рассматривая отдельные вычисления.
Инфраструктура для полиэдральной компиляции выполняет построение и оптимизацию полиэдрального представления программы, используя промежуточное представление программного кода, генерируемого компилятором. Промежуточное представление создается фронтендом компилятора с целью дальнейшей оптимизации и генерации кода для выбранной архитектуры [15].
В работах [9, 10] авторами показано, что аналитическое моделирование МММ, выполняемое BLIS, и автоматизированная оптимизация и распараллеливание, выполняемые над полиэдральным представлением, могут быть использованы для достижения производительности реализаций, созданных для случая свертки тензоров ТС. Отметим, что MVM является важным частным случаем ТС.
На рис. 1 изображен комплекс автоматизированной оптимизации обобщенных тензорных операций (АООТО), представленный в работах [9, 10]. Комплекс АООТО имеет три основные части: фронтенд, выполняющий компиляцию программы и построение промежуточного представления; инфраструктуру для построения полиэдрального представления и его модификаций, например, для распознавания ТС и ее оптимизации (на рис. 1 выделено цветом); бэкенд, выполняющий оптимизацию промежуточного кода программы, а также генерацию ассемблерного кода для целевой архитектуры. Такая структура позволяет обеспечить тестируемость описанных частей и их взаимозаменяемость.
В качестве фронтенда, выполняющего компиляцию в комплексе АООТО, использован фронтенд Clang [16]. Комплекс оптимизаций Polly [15] применен в качестве инфраструктуры для построения полиэдрального представления, а также в качестве основы для реализации распознавания ТС и его оптимизации. Внешняя библиотека LLVM Core [16] использована в качестве бэкенда, выполняющего оптимизацию промежуточного представления программы, а также генерацию ассемблерного кода для целевой архитектуры.
Существует множество подходов для оптимизации ТС по времени выполнения. В частности, реализация ТС может быть оптимизирована, используя последовательности оптимизаций циклов [17]. Для оптимизации ТС могут быть использованы высокопроизводительные готовые реализации. Фреймворки GETT (The general matrix multiplication (GEMM)-like
Tensor-Tensor multiplication) [18] и TBLIS (Tensor-Based Library Instantiation Software) [19] используют микро-ядра МММ. Микро-ядра являются частями готовой реализации МММ, содержащие векторные инструкции. Фреймворк Polly использует микро-ядра МММ и генерирует их автоматизировано.
Программа
Инфраструктура для полиэдральных оптимизаций
Построение полиэдрального пре д с тав ле ния
N Определение
Распознавание оптимальных
ТС значений
V ) параметров ТС
( S
-> Оптимизация ТС
)
Фронте нд
Построение про ме жуто чно го представления
Промежуточное
представление
Оптимизация кода
Генерация кода для целевой архитектуры
Ассемблерный код
Рис. 1. Комплекс программ
Polly не поддерживает оптимизацию MVM. В ближайшее время предполагается внедрение описанного алгоритма ААО MVM в Polly. Новые инструкции полиэдрального представления могут быть сгенерированы для выполнения предвыборки данных в кэш первого уровня выбранных итераций циклов. Для нетранспонированных матриц в MVM нужно сгенерировать и векторизовать редукцию элементов вычисляемого вектора [5, 20].
2. Численные эксперименты
2.1. Постановка обратной задачи гравиметрии и метод решения
Рассмотрим трехмерную структурную обратную задачу гравиметрии о восстановлении поверхности раздела между средами.
Задача состоит в нахождении функции ( = ((х,у), описывающей поверхность раздела сред постоянной плотности о\ и 02 по гравитационному полю Ад(х,у,0), измеренному на некоторой площади земной поверхности, и скачку плотности на границе раздела Асг = ci — (72- Пусть z = h — асимптотическая плоскость для данной поверхности раздела (
такая, что lim С(ж,у) = h. х,у—>-±оо
Задача описывается двумерным нелинейным интегральным уравнением первого рода [21]:
1
JL±0
— ОО —ОО
\/{{х ~ ж')2 + {у- УО2 + С2(У, У'))
V((х - х')2 + (у - у')2 + h2)
>dx'dy' = Ад(х,у,0),
(1)
где / — гравитационная постоянная.
Предположим, что значения поля вне некоторой прямоугольной области П = {(х,у) : a^x^b,c^y^d} близки к нулю. Проведем дискретизацию области П на сетке N х М : {xiiVj)i i = l,...,N, j = 1,..., М с шагами Аж и Ау.
После аппроксимации интегрального оператора по квадратурным формулам получим уравнение
/АаАхАух
N М Г
*ЕЕ
\Дх1
1
-xu)2+(yj-yv)2+C2(xi, yj)
1
л/(Xi - Хи)2 + (yj - yv)2 + К2
Ag(xu, yv, 0),
(la)
u = 1,M, v = 1,N.
Правую часть Ag(xu,yv,0) представим в виде вектора
F = {F(v-!)M+u} = {Ag(xu,yv, 0)}, и = 1,М, v = 1,N
и искомую поверхность представим в виде вектора
2 = {z(j-i)M+i} = {C(xi,yj), i = 1,М, j = 1,N}.
Тогда систему (la) запишем в операторном виде
Mz) = F. (2)
Для решения системы (2) будем использовать метод Левенберга—Марквардта [22]
zk+1 = zk-7
А\гк)У A\zk)+al
-i
xA'(zky(A'(zk)-F),
(3)
где zk — приближенное решение на к-ой итерации, A,(zk) — производная Фреше оператора А в точке ик, I — единичный оператор, 7 — демпфирующий множитель, а — параметр регул яриз ации.
В качестве начального приближения используется z° = h. Критерием останова является условие ||A(,2fc) — КЦ/ЦКЦ ^ £\ при некотором £\ < 0 .
Вместо выполнения трудоемкой процедуры обращения матрицы воспользуемся следующим приемом. На каждом шаге итерационного метода Левенберга—Марквардта [22] решаем систему линейных алгебраических уравнений
Bz = Ъ, (4)
где
В= (Л'(/))Т A!(zk) + cd, b= \(A!{zk)Y A!{zk) + al] zk - 'yA'(zk)J(A'(zk) - F).
Для решения (4) используется метод минимальных невязок [22]
/+1_ / (B(Bzl-b),Bzl-b)fTyJ
z = z —
II B(Bz1-.
-{В z1-b),
(5)
где z1 — приближенное решение на I-ой итерации метода минимальных невязок.
В качестве начального приближения используется z° = 0. Критерием останова является условие ||Bzl — б||/||Ь|| ^ £2 при некотором £2 < 0 .
Замечание. Метод Левенберга—Марквардта (3) и метод минимальных невязок (5) на каждом шаге итерационного процесса (3) требуют выполнения значительного количества операций MVM. Поэтому предложенный алгоритм автоматической оптимизации MVM может быть использован для автоматизированной и эффективной оптимизации соответствующей программы.
2.2. Результаты экспериментов
В данном разделе проведено сравнение производительности нескольких реализаций программного кода для решения обратной задачи гравиметрии (раздел 2.1) с использованием библиотек линейной алгебры Intel MKL, BLIS, OpenBLAS с реализацией, полученной с использованием предлагаемого алгоритма ААО MVM (см. раздел 1).
Для выполнения экспериментов использовались два 18-ядерных процессора Intel Хеоп Е5-2697 v4 суперкомпьютера «Уран». Для распараллеливания программы применялись директивы стандарта ОрепМР и группы из 36 потоков.
Intel MKL --Q-- OpenBLAS
—0— BLIS 0 Model
Размер квадратной сетки
Рис. 2. Ускорение, полученное для обратной задачи гравиметрии
Определим ускорение S = Tseriai/Tparauei, где Tseriai - время выполнения последовательного кода программы с применением стандартных оптимизаций компилятора третьего уровня (03), Tparaiiei - время выполнения параллельного кода программы с использованием библиотечных процедур либо кода, оптимизированного с помощью предлагаемого алгоритма.
На рис. 2 для 36 ядер процессора Intel Хеоп Е5-2697 v4 суперкомпьютера «Уран» изображены коэффициенты ускорения для решения обратной задачи гравиметрии методом Левенберга—Марквардта на квадратных сетках с размерами, изменяющимися от 32 до
160 с шагом 4. Представленные коэффициенты ускорения получены для разработанного алгоритма А АО MVM, а также для библиотек BLIS, OpenBLAS и Intel MKL.
В большинстве случаев наилучшие результаты получены с использованием библиотеки Intel MKL. Результаты алгоритма ААО MVM сравнимы с библиотечными для всех рассмотренных размеров сеток. В среднем они отличаются не более чем на 1% и превосходят результаты библиотек OpenBLAS и BLIS.
На основании результатов численных экспериментов можно сделать вывод, что реализация программы, использующей описанный метод оптимизации MVM, сравнима по производительности с реализациями, использующими библиотеки BLIS, OpenBLAS и Intel MKL, оптимизированные вручную. Результаты численных экспериментов показывают, что значения параметров реализации, полученные с использованием аналитического моделирования (раздел 1), не зависят от размерности задачи.
Заключение
В данной работе представлен алгоритм автоматизированной оптимизации матричновекторного произведения по времени выполнения, использующийся на этапе компиляции. Алгоритм ААО MVM основан на оптимизации полиэдрального представления программы и моделировании вычислений на гипотетическом многоядерном процессоре. В отличие от подходов, основанных на ручной настройке и автонастройке, описанный алгоритм может применяться для создания новых реализаций матрично-векторного произведения в условиях недоступности целевой архитектуры процессора и ограниченности времени выполнения.
Алгоритм использован для оптимизации программного кода, реализующего решение структурной обратной задачи гравиметрии о нахождении поверхности раздела сред методом Левенберга—Марквардта. Проведено сравнение производительности полученной реализации с реализациями на основе оптимизированных библиотек линейной алгебры Intel MKL, BLIS, OpenBLAS. Результаты численных экспериментов показывают сравнимость предложенного алгоритма по эффективности с подходами, созданными с использованием эмпирического поиска при доступе к целевым архитектурам процессоров. В дальнейшем планируется внедрение описанного алгоритма в комплекс оптимизаций Polly с целью использования в процессе компиляции, выполняемой фронтендом Clang.
Перспективным направлением будущих исследований является применение аналитического моделирования для оптимизации матрично-векторного произведения матриц специального вида (матрица Теплица и др.). В основе таких подходов может использоваться симметрия элементов и другие свойства матриц.
Литература
1. Yu J., Lukefahr A., Palframan D., et al. Scalpel: Customizing DNN Pruning to the Underlying Hardware Parallelism // SIGARCH Computer Architecture News. 2017. Vol. 45, no. 2, P. 548-560. DOI: 10.1145/3140659.3080215.
2. Yang X., Parthasarathy S., Sadayappan P. Fast Sparse Matrix-vector Multiplication on GPUs: Implications for Graph Mining // Proceedings of the VLDB Endowment. 2011. Vol. 4, no. 4. P. 231-242. DOI: 10.14778/1938545.1938548.
3. Kaushik D., Gropp W., Minkoff M., et al. Improving the Performance of Tensor Matrix Vector Multiplication in Cumulative Reaction Probability Based Quantum Chemistry Codes // High Performance Computing (HiPC 2008). Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
Р. 120-130. DOI: 10.2172/928654.
4. Martyshko P.S., Akimova E.N., Misilov V.E. Solving the structural inverse gravity problem by the modified gradient methods // Izvestiya, Physics of the Solid Earth. 2016. Vol. 52, no. 5. P. 704-708. DOI: 10.1134/S1069351316050098.
5. Hassan S.A., Mahmoud M.M., Hemeida A., et al. Effective Implementation of Matrix-Vector Multiplication on Intel’s AVX Multicore Processor // Computer Languages, Systems & Structures. 2018. Vol. 51, P. 158-175. DOI: 10.1016/j.cl.2017.06.003.
6. Liang J., Zhang Y. Optimization of GEMV on Intel AVX processor // International Journal of Database Theory and Application. 2016. Vol. 9. P. 47-60. DOI: 10.14257/ijdta.2016.9.2.06.
7. Low T.M., Igual F.D., Smith T.M. Quintana-Orti E.S., Analytical Modeling Is Enough for High-Performance BLIS // ACM Transactions on Mathematical Software. 2016. Vol. 43, no. 2. P. 12:1-12:18. DOI: 10.1145/2925987.
8. Frison G. Algorithms and Methods for High-Performance Model Predictive Control / / Technical University of Denmark. 2016. 345 c.
9. Akimova E.N., Gareev R.A. Algorithm of Automatic Parallelization of Generalized Matrix Multiplication // CEUR Workshop Proceedings. 2017. Vol. 2005. P. 1-10.
10. Gareev R., Grosser T., Kruse M. High-Performance Generalized Tensor Operations: A Compiler-Oriented Approach // ACM Transactions on Architecture and Code Optimization. 2018. Vol. 15, no. 3. P. 34:1-34:27. DOI: 10.1145/3235029.
11. Intel. Intel Math Kernel Library (Intel MKL). URL: https: //software. intel. com/en-us/ mkl (дата обращения: 27.10.2019).
12. Van Zee F.G., Van De Geijn R.A. BLIS: A Framework for Rapidly Instantiating BLAS Functionality // ACM Transactions on Mathematical Software. 2015. Vol. 41, no. 3. P. 14:1-14:33. DOI: 10.1145/2764454.
13. Xianyi Z., Qian W., Yunquan Z. Model-driven level 3 BLAS performance optimization on Loongson ЗА processor // Parallel and Distributed Systems (ICPADS), 2012 IEEE 18th International Conference on Parallel and Distributed Systems, IEEE. 2012. P. 684-691. DOI: 10.1109/icpads.2012.97.
14. Feautrier P., Lengauer C. Polyhedron Model // Encyclopedia of Parallel Computing. Springer US, Boston, MA. 2011. P. 1581-1592. DOI: 10.1007/978-0-387-09766-4_502.
15. Grosser T., Grófflinger A., Lengauer C. Polly—Performing polyhedral optimizations on a low-level intermediate representation // Parallel Process. Lett. 2012. Vol. 22. no. 4. DOI: 10.1142/S0129626412500107.
16. Lattner C. LLVM: An Infrastructure for Multi-Stage Optimization // Master’s Thesis. 2002. URL: http://llvm.cs.uiuc.edu (дата обращения: 27.10.2019).
17. Apra E., Klemm M., Kowalski K. Efficient Implementation of Many-Body Quantum Chemical Methods on the Intel@Xeon Phi™Coprocessor // International Conference for High Performance Computing, Networking, Storage and Analysis. 2014. P. 674-684. DOI: 10.1109/SC.2014.60
18. Springer P., Bientinesi P. Design of a high-performance GEMM-like tensor-tensor multiplication / / ACM Trans. Math. Softw. 2018. Vol. 44, no. 3, Article 28. DOI: 10.1145/3157733.
Scientific Computing. 2016. Vol. 40. DOI: 10.1137/16M108968X
20. Doerfert J., Streit К., Hack S., et al. Polly’s polyhedral scheduling in the presence of reductions // arXiv preprint. 2015. arXiv: 1505.07716.
21. Numerov B.V. Interpretation of gravitational observations in the case of one contact surface // Doki. Akad. Nauk SSSR. 1930. P. 569-574.
22. Vasin V.V., Eremin 1.1. Operators and iterative processes of Fejer type: theory and
applications // Walter de Gruyter. 2009. Vol. 53. 155 p. DOI: 10.1515/9783110218190.
Акимова Елена Николаевна, д.ф.-м.н., доцент, ведущий научный сотрудник Института математики и механики им. Н.Н. Красовского Уральского отделения РАН (Екатеринбург, Российская Федерация); профессор кафедры информационных технологий и систем управления, Институт радиоэлектроники и информационных технологий-РтФ, УрФУ имени первого Президента России Б.Н. Ельцина (Екатеринбург, Российская Федерация)
Гареев Роман Альбертович, аспирант, Институт радиоэлектроники и информационных технологий-РтФ, УрФУ имени первого Президента России Б.Н. Ельцина (Екатеринбург, Российская Федерация)
DOI: 10.14529/cmse200105
APPLICATION OF ANALYTICAL MODELING OF MATRIX-VECTOR MULTIPLICATION ON MULTICORE PROCESSORS
© 2020 E.N. Akimova1,2, R.A. Gareev1
Many areas of study and their applications such as machine learning, data mining, quantum chemistry, mathematical physics, and high-performance computing require effective implementation of matrix-vector multiplication. In this work, we present an overview of the algorithm for automatic optimization of matrix-vector multiplication. The algorithm models computations on the hypothetical multicore processor, which is introduced by the authors, and applies polyhedral modeling. In comparison with methods, which rely on manual tuning and auto-tuning, the algorithm can be utilized when the execution time is a critical factor and the target platform is not accessible. We apply the approach to optimize an implementation of the solution of the inverse gravimetry problem of finding an interface between the layers, which uses the iterative Levenberg—Marquardt method. The performance of the obtained implementation is compared with the performances produced by implementations, which are based on MKL, BLIS, and OpenBLAS. Results of the experimental evaluation show that the considered approach is comparable with the approaches, which are created on target architectures using manual tuning.
Keywords: compilers, linear algebra, matrix-vector operations, analytical modeling, inverse gravimetry problem.
1 Ural Federal University (Mira 19, Yekaterinburg, 620002 Russia),
2N.N. Krasovskii Institute of Mathematic sand Mechanics, UrB RAS (S. Kovalevskaya 16, Yekaterinburg, 620990 Russia)
E-mail: aenl5@yandex.ru, gareevroman@gmail.com Received: 31.01.2020
FOR CITATION
Akimova E.N., Gareev R.A. Application of Analytical Modeling of Matrix-Vector Multiplication on Multicore Processors. Bulletin of the South Ural State University. Series: Computational Mathematics and Software Engineering. 2020. Vol. 9, no. 1. P. 69-82. (in Russian) DOI: 10.14529/cmse200105.
This paper is distributed under the terms of the Creative Commons Attribution-Non Commercial 3.0 License which permits non-commercial use, reproduction and distribution of the work without further permission provided the original work is properly cited.
References
1. Yu J., Lukefahr A., Palframan D., et al. Scalpel: Customizing DNN Pruning to the Underlying Hardware Parallelism. SIGARCH Computer Architecture News. 2017. Vol. 45, no. 2. P. 548-560. DOI: 10.1145/3140659.3080215.
2. Yang X., Parthasarathy S., Sadayappan P. Fast Sparse Matrix-vector Multiplication on GPUs: Implications for Graph Mining. Proceedings of the VLDB Endowment. 2011. Vol. 4, no. 4. P. 231-242. DOI: 10.14778/1938545.1938548.
3. Kaushik D., Gropp W., Minkoff M., et al. Improving the Performance of Tensor Matrix Vector Multiplication in Cumulative Reaction Probability Based Quantum Chemistry Codes. High Performance Computing (HiPC 2008). Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. P. 120-130. DOI: 10.2172/928654.
4. Martyshko P.S., Akimova E.N., Misilov V.E. Solving the structural inverse gravity problem by the modified gradient methods. Izvestiya, Physics of the Solid Earth. 2016. Vol. 52, no. 5. P. 704-708. DOI: 10.1134/S1069351316050098.
5. Hassan S.A., Mahmoud M.M., Hemeida A., et al. Effective Implementation of Matrix-Vector Multiplication on Intel’s AVX Multicore Processor. Computer Languages, Systems & Structures. 2018. Vol. 51. P. 158-175. DOI: 10.1016/j.cl.2017.06.003.
6. Liang J., Zhang Y. Optimization of GEMV on Intel AVX processor. International Journal of Database Theory and Application. 2016. Vol. 9. P. 47-60. DOI: 10.14257/ijdta.2016.9.2.06.
7. Low T.M., Igual F.D., Smith T.M. Quintana-Orti E.S., Analytical Modeling Is Enough for High-Performance BLIS. ACM Transactions on Mathematical Software. 2016. Vol. 43, no. 2. P. 12:1-12:18. DOI: 10.1145/2925987.
8. Frison G. Algorithms and Methods for High-Performance Model Predictive Control. Technical University of Denmark. 2016. 345 p.
9. Akimova E.N., Gareev R.A. Algorithm of Automatic Parallelization of Generalized Matrix Multiplication. CEUR Workshop Proceedings. 2017. Vol. 2005. P. 1-10.
10. Gareev R., Grosser T., Kruse M. High-Performance Generalized Tensor Operations: A Compiler-Oriented Approach. ACM Transactions on Architecture and Code Optimization. 2018. Vol. 15, no. 3. P. 34:1-34:27. DOI: 10.1145/3235029.
11. Intel. Intel Math Kernel Library (Intel MKL). Available at: https: //software. intel. com/ en-us/mkl (accessed: 27.10.2019).
12. Van Zee F.G., Van De Geijn R.A. BLIS: A Framework for Rapidly Instantiating BLAS Functionality. ACM Transactions on Mathematical Software. 2015. Vol. 41, no. 3. P. 14:1-14:33. DOI: 10.1145/2764454.
13. Xianyi Z., Qian W., Yunquan Z. Model-driven level 3 BLAS performance optimization on Loongson ЗА processor. Parallel and Distributed Systems (ICPADS), 2012 IEEE 18th International Conference on Parallel and Distributed Systems, IEEE. 2012. P. 684-691. DOI: 10.1109/icpads.2012.97.
14. Feautrier P., Lengauer C. Polyhedron Model. Encyclopedia of Parallel Computing. Springer US, Boston, MA. 2011. P. 1581-1592. DOI: 10.1007/978-0 387-09766-4_502.
15. Grosser T., Grofflinger A., Lengauer C. Polly—Performing polyhedral optimizations on a low-level intermediate representation. Parallel Process. Lett. 2012. Vol. 22, no. 4. DOI: 10.1142/S0129626412500107.
16. Lattner C. LLVM: An Infrastructure for Multi-Stage Optimization. Master’s Thesis. 2002. Available at: http://llvm.cs.uiuc.edu (accessed: 27.10.2019).
17. Apra E., Klemm M., Kowalski K. Efficient Implementation of Many-Body Quantum Chemical Methods on the Intel@Xeon Phi™Coprocessor. International Conference for High Performance Computing, Networking, Storage and Analysis. 2014. P. 674-684. DOI: 10.1109/SC.2014.60
18. Springer P., Bientinesi P. Design of a high-performance GEMM-like tensor-tensor multiplication. ACM Trans. Math. Softw. 2018. Vol. 44, no. 3, Article 28. DOI: 10.1145/3157733.
19. Matthews D. High-performance tensor contraction without BLAS. SIAM Journal on Scientific Computing. 2016. Vol. 40. DOI: 10.1137/16M108968X.
20. Doerfert J., Streit K., Hack S., et al. Polly’s polyhedral scheduling in the presence of reductions. arXiv preprint. 2015. arXiv:1505.07716.
21. Numerov B.V. Interpretation of gravitational observations in the case of one contact surface. Doki. Akad. Nauk SSSR. 1930. P. 569-574.
22. Vasin V.V., Eremin I.I. Operators and iterative processes of Fejer type: theory and
applications. Walter de Gruyter. 2009. Vol. 53. 155 p. DOI: 10.1515/9783110218190.
