Перемножение плотных неособенных матриц (последовательный вещественный вариант)

Основные авторы описания: А.В.Фролов, Вад.В.Воеводин (раздел 2.2), А.М.Теплов (раздел 2.4)

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Перемножение матриц - одна из базовых задач в алгоритмах линейной алгебры, широко применяется в большом количестве разных методов. Здесь мы рассмотрим умножение [math]C = AB[/math]  плотных неособенных матриц (последовательный вещественный вариант), то есть тот вариант, где никак не используются ни специальный вид матрицы, ни ассоциативные свойства операции сложения.

1.2 Математическое описание алгоритма

Исходные данные: плотная матрица [math]A[/math] (элементы [math]a_{ij}[/math]), плотная матрица [math]B[/math] (элементы [math]b_{ij}[/math]).

Вычисляемые данные: плотная матрица [math]C[/math] (элементы [math]c_{ij}[/math]).

Формулы метода:

[math] \begin{align} c_{ij} = \sum_{k = 1}^{n} a_{ik} b_{kj}, \quad i \in [1, m], \quad i \in [1, l]. \end{align} [/math]

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

1.3 Вычислительное ядро алгоритма

Вычислительное ядро перемножения плотных неособенных матриц можно составить из множественных (всего их [math]l[/math]) вычислений умножения матрицы [math]A[/math] на столбцы матрицы [math]B[/math], или (при более детальном рассмотрении), из множественных (всего их [math]ml[/math]) скалярных произведений строк матрицы [math]A[/math] на столбцы матрицы [math]B[/math]:

[math]\sum_{k = 1}^{n} a_{ik} b_{kj}[/math]

в режиме накопления или без него, в зависимости от требований задачи.

1.4 Макроструктура алгоритма

Как уже записано в описании ядра алгоритма, основную часть умножения матриц составляют множественные (всего [math]ml[/math]) вычисления скалярных произведений строк матрицы [math]A[/math] на столбцы матрицы [math]B[/math]

[math]\sum_{k = 1}^{n} a_{ik} b_{kj}[/math]

в режиме накопления или без него.

1.5 Схема реализации последовательного алгоритма

Для всех [math]i[/math] от [math]1[/math] до [math]m[/math] и для всех [math]j[/math] от [math]1[/math] до [math]l[/math] выполняются

[math]c_{ij} = \sum_{k = 1}^{n} a_{ik} b_{kj}[/math]

Особо отметим, что вычисления сумм вида [math]\sum_{k = 1}^{n} a_{ik} b_{kj}[/math] производят в режиме накопления прибавлением к текущему (временному) значению вычисляемого элемента матрицы [math]c_{ij}[/math] произведений [math]a_{ik} b_{kj}[/math] для [math]k[/math] от [math]1[/math] до [math]n[/math], c возрастанием [math]k[/math], вначале все элементы инициализируются нулями. При суммировании "по убыванию" общая схема принципиально не отличается и потому нами не рассматривается. Другие порядки выполнения суммирования приводят к изменению параллельных свойств алгоритма и будут рассматриваться нами в отдельных описаниях.

1.6 Последовательная сложность алгоритма

Для умножения двух квадратных матриц порядка [math]n[/math] (т.е. при [math]m=n=l[/math]) в последовательном (наиболее быстром) варианте требуется:

• по [math]n^3[/math] умножений и сложений.

Для умножения матрицы размером [math]m[/math] строк на [math]n[/math] столбцов на матрицу размером [math]m[/math] строк на [math]n[/math] столбцов в последовательном (наиболее быстром) варианте требуется:

• по [math]mnl[/math] умножений и сложений.

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

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

1.7 Информационный граф

Опишем граф алгоритма как аналитически, так и в виде рисунка.

Граф алгоритма умножения плотных матриц состоит из одной группы вершин, расположенной в целочисленных узлах трёхмерной области, соответствующая ей операция [math]a+bc[/math].

Естественно введённые координаты области таковы:

• [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]m[/math], принимая все целочисленные значения;
• [math]j[/math] — меняется в диапазоне от [math]1[/math] до [math]l[/math], принимая все целочисленные значения;
• [math]k[/math] — меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения.

Аргументы операции следующие:

• [math]a[/math]:
  • при [math]k = 1[/math] константа [math]0[/math];
  • при [math]k \gt 1[/math] — результат срабатывания операции, соответствующей вершине с координатами [math]i, j, k-1[/math];
• [math]b[/math] — элемент входных данных, а именно [math]a_{ik}[/math];
• [math]c[/math] - элемент входных данных [math]b_{kj}[/math];

Результат срабатывания операции является:

• при [math]k \lt n[/math] - промежуточным данным алгоритма;
• при [math]k = n[/math] - выходным данным [math]c_{ij}[/math].

Рисунок 1. Умножение плотных матриц с отображением выходных данных

Интерактивное изображение графа алгоритма без входных и выходных данных для случая перемножения двух квадратных матриц порядка 3 и 4

1.8 Ресурс параллелизма алгоритма

Для алгоритма умножения квадратных матриц порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:

• по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]n^2[/math] операций).

Для умножения матрицы размером [math]m[/math] строк на [math]n[/math] столбцов на матрицу размером [math]n[/math] строк на [math]l[/math] столбцов в последовательном (наиболее быстром) варианте требуется:

• по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]ml[/math] операций).

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

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

1.9 Входные и выходные данные алгоритма

Входные данные: матрица [math]A[/math] (элементы [math]a_{ij}[/math]), матрица [math]B[/math] (элементы [math]b_{ij}[/math])).

Объём входных данных: [math]mn+nl[/math]

Выходные данные: матрица [math]C[/math] (элементы [math]c_{ij}[/math]).

Объём выходных данных: [math]ml[/math]

1.10 Свойства алгоритма

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

При этом вычислительная мощность алгоритма умножения матриц, как отношение числа операций к суммарному объему входных и выходных данных – линейно.

При этом алгоритм умножения матриц полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

В простейшем варианте алгоритм умножения матриц на Фортране можно записать так:

         
	DO  I = 1, M
        DO  J = 1, L
		S = 0.
		DO  K = 1, N
			S = S + DPROD(A(I,K), B(K,J))
		END DO	
	        C(I, J) = S
        END DO
	END DO

При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.

2.2 Локальность данных и вычислений

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

Рисунок 2. Профили обращений для 6 вариантов перемножения матриц

На рис.2 представлены 6 профилей обращений для различных вариантов классического перемножения матриц (в зависимости от выбранного порядка циклов). На каждом профиле хорошо выделяются фрагменты профилей для 3-х массивов, используемых в программе. При этом порядок обращений к разным массивам во всех 6 вариантах один и тот же; т.е. взаимодействие между массивами везде устроен одинаково. В этом случае отличия в локальности задаются внутренним строением фрагментов профиля для каждого массива в отдельности.

Исходя из исходного кода, можно увидеть, что всего встречается 6 разных видов фрагментов для отдельных массивов (эти виды выделены на рис.2 зеленым). Стоит сделать два уточнения: 1) профиль для результирующего массива С всегда в 2 раза больше, поскольку к элементам этого массива всегда происходят по два обращения подряд; 2) если во внутреннем цикле используется элементы массива, то обращение к этому элементу выносится за пределы внутреннего цикла (в цикле используется скалярная переменная вместо него).

Также заранее отметим, что рисунки по профилям получаются следующим образом: строится рисунок по общему профилю, после чего из него оставляется только часть для рассматриваемого массива. Доля обращений к одному массиву может меняться, поэтому и частота обращений на каждом рисунке может в значительной степени отличаться.

Рассмотрим подробнее каждый из данных 6 фрагментов на примере массива С.

Фрагмент 1. На рис.3 показано начало фрагмента 1 (здесь и далее начало фрагмента соответствует выделенной оранжевой области на рис.2). Можно увидеть, что данный фрагмент устроен достаточно просто: выполняется перебор некоторого блока элементов массива, затем данный перебор циклически повторяется. После этого происходит переход к следующему блоку элементов массива, и выполняется тот же перебор в цикле.

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

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

Рисунок 3. Начало фрагмента 1

Рисунок 4. Начало фрагмента 1, первые 3 итерации

Фрагмент 2. На рис.5 представлено начало фрагмента 2. Видно, что данный фрагмент также представляет собой перебор элементов массива в цикле, однако в данном случае на каждой итерации цикла перебираются сразу все элементы массива, а не только элементы одного блока. При более подробном рассмотрении итерации (рис.6) видно, что здесь также происходит именно последовательный перебор, и также к каждому элементу обращение происходит дважды (поскольку это массив С).

Рисунок 5. Начало фрагмента 2

Рисунок 6. Начало фрагмента 2, часть первой итерации

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

Фрагмент 3. Начало данного фрагмента представлено на рисунке 7. Данный профиль является последовательным перебором всех элементов массива, однако более подробное рассмотрение показывает, что к каждому элементу происходит обращение дважды, при этом между обращениями происходит достаточно много других обращений. Это связано с указанным ранее замечанием касательно вынесения повторяющегося обращения за предел внутреннего цикла.

Рисунок 7. Начало фрагмента 3

Рисунок 8. Начало фрагмента 3, первые несколько обращений

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

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

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

Рисунок 9. Начало фрагмента 4

Фрагмент 5. В отличие от предыдущих фрагментов, исходя из визуализации начала данного фрагмента (рис.9), сложно сделать выводы относительно структуры обращения в память. Однако более подробное рассмотрение (рис.10) показывает, что данный профиль практически идентичен профилю предыдущего фрагмента. Более того, фрагмент 5 на самом деле состоит из итерационно повторяющихся фрагментов 4.

Рисунок 10. Начало фрагмента 5

Рисунок 11. Начало фрагмента 5, первые несколько итераций

Данный фрагмент обладает по сравнению с фрагментом 4 практически такой же пространственной локальностью, однако более высокой временной локальностью. Причина в обоих случаях одна – в данном случае тот же набор обращений повторяется несколько раз.

Фрагмент 6. Данный фрагмент, представленный на рис.12 и рис.13, так же сильно напоминает фрагмент 5, но с одним отличием. Во фрагменте 5 несколько раз повторяется фрагмент 4, в рамках которого выполняется несколько итераций, и на каждой итерации внутри фрагмента 4 используются разные элементы (поскольку первый элемент сдвигается на 1). В случае данного фрагмента выполняются все те же итерации, но в перемешанном виде. Сначала выполняются все итерации, начинающиеся с одного и того же элемента; затем все итерации, который начинаются с элемента, чей индекс идет следующим; и т.д.

Такой фрагмент обладает более высокой и пространственной, и временно́й локальностью по сравнению с фрагментом 5, поскольку обращения в одним и тем же элементам расположены ближе друг к другу в профиле.

Рисунок 12. Начало фрагмента 6

Рисунок 13. Начало фрагмента 6, первые несколько итераций

В итоге можно сказать, что наиболее низкой локальностью в целом обладают фрагменты 5 и 6. Фрагмент 4 также обладает низкой локальностью, однако содержит значительно меньшее число обращений, а значит, привносит меньший вклад в общий профиль обращений. Это позволяет определить, что наименее эффективными с точки зрения работы с памятью являются варианты kji и jki, поскольку каждый из них содержит фрагменты 5 и 6. Далее идут варианты ijk и jik – в них по одному такому фрагменту. Самые лучшие варианты – ikj и kij.

2.2.1.2 Количественная оценка локальности

Основные фрагменты реализации, на основе которых были получены количественные оценки, приведены здесь (функция Kernel***, где *** - порядок циклов (напр., KernelIJK) ). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

На рис.14 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что, как и предполагалось наиболее эффективными являются варианты kij,ikj. Заметно хуже результат у вариантов jik,ijk, при этом эти варианты практически равны между собой. Самый плохой результат ожидаемо у вариантов jki,kji.

Рисунок 14. Сравнение значений оценки daps

Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.

На рис.15 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Относительное сравнение результатов по cvg в целом повторяет результаты daps.

Рисунок 15. Сравнение значений оценки cvg

2.4 Масштабируемость алгоритма и его реализации

2.4.2 Масштабируемость реализации алгоритма

Рисунок 16. Параллельная реализация произведения матриц Максимальная производительность.

Рисунок 17. Параллельная реализация произведения матриц Максимальная эффективность.

Набор изменяемых параметров запуска реализации алгоритма и границы значений параметров алгоритма:

• число процессоров [4 : 1024]
• размерность матрицы [1024 : 20480]

Эффективность выполнения реализации алгоритма

• Минимальная эффективность 4,71%
• Максимальная эффективность 31,72%

Оценка масштабируемости

• По числу процессов: -0.0436 – при увеличении числа процессов эффективность убывает достаточно интенсивно на всей рассмотренной области изменений параметров запуска. Уменьшение эффективности на рассмотренной области работы параллельной программы звязано с увеличением числа пересылок с ростом числа процессов и как следствие ростом накладных расходов на организацию вычислений. Присутствует область, на которой при увеличении числа процессов эффективность возрастает, но при дальнейшем росте продолжает снижаться. Это объясняется декомпозицей данных, при которой наступает момент, когда размер матрицы позволяет блокам укладываться в КЭШ-память. Так же это подтверждает проявление этого явления, но со смещением по числу процессов, и при увеличении вычислительной сложности задачи.
• По размеру задачи: -0.0255 – при увеличении размера задачи эффективность в целом уменьшается по рассматриваемой области, хотя и менее интенсивно, чем при увеличении числа процессов. Снижение эффективности объясняется тем, что при росте вычислительной сложности существенно возрастают объемы передаваемых данных. Присутствует область возрастания эффективности, на всех рассмотренных размерах матрицы. Это объясняется тем, что при малом размере задачи данные хорошо укладываются в КЭШ память, что приводит к высокой эффективности работы приложения при малом размере задачи. При дальнейшем увеличении размера эффективность уменьшается при выходе за границы КЭШ-памяти.
• По двум направлениям: -0.000968 – при рассмотрении увеличения, как вычислительной сложности, так и числа процессов по всей рассмотренной области значений уменьшается, интенсивность уменьшения эффективности не очень высока. В совокупности с тем фактом, что разница между максимальной и минимальной эффективностью на рассмотренной области значений параметров составляет почти 25% говорит о том, что уменьшение эффективности по всей области довольно равномерное, но интенсивно лишь в не очень больших участках по площади. На остальной области значений параметров изменения эффективности не столь значительны и находятся на приблизительно одном и том же уровне.

Реализация алгоритма на языке C.