Участник:Diana Pimenova/Алгоритм Фокса умножения матриц

Автор: Д.В.Пименова

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

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

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

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

Исходные данные: квадратная матрица [math]\displaystyle{ A }[/math] (состоит из блоков [math]\displaystyle{ A_{ij} }[/math]), квадратная матрица [math]\displaystyle{ B }[/math] (состоит из блоков [math]\displaystyle{ B_{ij} }[/math]).

Вычисляемые данные: квадратная матрица [math]\displaystyle{ C }[/math] (состоит из блоков [math]\displaystyle{ C_{ij} }[/math]).

При этом каждый блок [math]\displaystyle{ C_{ij} }[/math] определяется как произведение соответствующих блоков матриц [math]\displaystyle{ A }[/math] и [math]\displaystyle{ B }[/math] :

[math]\displaystyle{ \begin{align} C_{ij} = \sum_{s = 0}^{q-1} A_{is} B_{sj}, \quad i \in [0, q-1],\quad j \in [0, q-1]. \end{align} }[/math]

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

Вычислительное ядро перемножения квадратных матриц методом Фокса можно составить из процедур вычисления всех элементов одного из блоков результирующей матрицы, т.е. из процедур множественных вычислений умножения блоков матрицы [math]\displaystyle{ A }[/math] на блоки матрицы [math]\displaystyle{ B }[/math]

[math]\displaystyle{ \begin{align} \sum_{s = 0}^{q-1} A_{is} B_{sj}, \quad i \in [0, q-1],\quad j \in [0, q-1]. \end{align} }[/math]

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

Как уже записано в описании ядра алгоритма, основную часть умножения матриц составляют множественные вычисления произведения блоков матрицы [math]\displaystyle{ A }[/math] на блоки матрицы [math]\displaystyle{ B }[/math]. Перемножать блоки будем стандартным методом, используя определение произведения матриц:

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

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

Разобьем исходные матрицы на блоки [math]\displaystyle{ A_{ij} }[/math],[math]\displaystyle{ B_{ij} }[/math] и [math]\displaystyle{ C_{ij} }[/math] соответственно. Обнулим блоки матрицы [math]\displaystyle{ C }[/math], предназначенной для хранения соответствующего результирующего произведения [math]\displaystyle{ AB }[/math]. Затем запускается цикл, в ходе которого выполняется последовательное умножение блоков матриц операндов и суммирование полученных результатов.

[math]\displaystyle{ \begin{align} С_{ij} = \sum_{s = 0}^{q-1} A_{is} B_{sj}, \quad i \in [0, q-1],\quad j \in [0, q-1]. \end{align} }[/math]

После завершения цикла мы получим матрицу с блоками [math]\displaystyle{ C_{ij} }[/math], соответствующую произведению [math]\displaystyle{ AB }[/math].

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

Рассмотрим квадратные матрицы размером [math]\displaystyle{ n \times n }[/math] , разбитые на блоки размера [math]\displaystyle{ \frac{n}{q} }[/math] . Итого у нас получается [math]\displaystyle{ q \times q }[/math]  блоков.

Для умножение данных матриц в последовательном варианте требуется по [math]\displaystyle{ n^3 }[/math]  умножений и сложений.

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

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

Сначала опишем идею параллельного умножения матриц методом Фокса:

• Каждой подзадаче [math]\displaystyle{ (i,j) }[/math] передаются блоки [math]\displaystyle{ A_{ij} }[/math], [math]\displaystyle{ B_{ij} }[/math] и обнуляются блоки [math]\displaystyle{ C_{ij} }[/math] на всех подзадачах.
• Для каждой строки [math]\displaystyle{ i }[/math] блок [math]\displaystyle{ A_{ij} }[/math] подзадачи [math]\displaystyle{ (i,j) }[/math] пересылается на все подзадачи той же строки [math]\displaystyle{ i }[/math] решетки.
• Полученные в результаты пересылок блоки [math]\displaystyle{ A'_{ij} }[/math], [math]\displaystyle{ B'_{ij} }[/math] каждой подзадачи [math]\displaystyle{ (i,j) }[/math] перемножаются и прибавляются к блоку [math]\displaystyle{ C_{ij} }[/math][math]\displaystyle{ : C_{ij} = C_{ij} + A'_{ij}B'_{ij} }[/math].
• блоки [math]\displaystyle{ B'_{ij} }[/math] каждой подзадачи [math]\displaystyle{ (i,j) }[/math] пересылаются подзадачам, являющимися соседями сверху в столбцах решетки подзадач (блоки подзадач из первой строки решетки пересылаются подзадачам последней строки решетки).

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

Определим вычислительную сложность данного алгоритма. Пусть все матрицы являются квадратными размера [math]\displaystyle{ n \times n }[/math] , количество блоков по горизонтали и вертикали являются одинаковым и равным [math]\displaystyle{ q }[/math]  (т.е. размер всех блоков равен [math]\displaystyle{ k \times k }[/math] , где [math]\displaystyle{ k = \frac{n}{q} }[/math] ), процессоры образуют квадратную решетку и их количество равно [math]\displaystyle{ p = q^{2} }[/math] .

Алгоритм Фокса требует для своего выполнения [math]\displaystyle{ q }[/math]  итераций, в ходе которых каждый процессор перемножает свои текущие блоки матриц [math]\displaystyle{ A }[/math]  и [math]\displaystyle{ B }[/math]  и прибавляет результаты умножения к текущему значению блока матрицы [math]\displaystyle{ C }[/math] . С учетом выдвинутых предположений общее количество выполняемых при этом операций будет иметь порядок [math]\displaystyle{ \frac{n^{3}}{p} }[/math] .

Определим количество вычислительных операций. Сложность выполнения скалярного умножения строки блока матрицы [math]\displaystyle{ A }[/math] на столбец блока матрицы [math]\displaystyle{ B }[/math] можно оценить как[math]\displaystyle{ 2\frac{n}{q} - 1 }[/math] . Количество строк и столбцов в блоках равно [math]\displaystyle{ \frac{n}{q} }[/math] и, как результат, трудоемкость операции блочного умножения оказывается равной [math]\displaystyle{ \frac{\frac{n^2}{p}}{\frac{2n}{q} - 1} }[/math]. Для сложения блоков требуется[math]\displaystyle{ \frac{n^{2}}{p} }[/math] операций.

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

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

Объём входных данных: [math]\displaystyle{ 2n^{2} }[/math]

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

Объём выходных данных: [math]\displaystyle{ n^{2} }[/math]

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

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

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

Возьмем квадратные матрицы размерностью [math]\displaystyle{ n \times n }[/math], где [math]\displaystyle{ 500\leqslant n \leqslant 2000 }[/math], а число процессов будем изменять от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ 128 }[/math]. В качестве примера приведем время выполнения программы в зависимости от числа процессов при фиксированном размере матриц [math]\displaystyle{ n = 2000 }[/math].

Число процессов	Время (с)
128	0.510
64	1.038
32	2.556
16	5.121
8	20.844
4	20.952
2	92.919

Далее рассмотрим зависимость времени выполнения программы от размера матриц. Число процессов также зафиксируем: возьмем [math]\displaystyle{ 64 }[/math] процесса.

Размер матриц	Время (с)
2000	1.038
1500	0.409
1000	0.115
750	0.049
500	0.015

Из приведенных результатов исследования видно, что система хорошо масштабируема.

Характеристики программно-аппаратной среды

Все вычисления были произведены на суперкомпьютере "Ломоносов".

Для компиляции был использован компилятор языка C++ GNU 4.4.6.

Вычисления производились в очереди test. Ограничений на лимит времени на узел наложено не было.

Существующие реализации алгоритма

Ссылки на некоторые существующие реализации этого алгоритма:

[1] - центр суперкомпьютерных технологий Нижегородского Государственного Университета.

[2] - ИНТУИТ