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

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

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

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

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

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

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

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

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

[math] \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]

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

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

[math] \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]

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

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

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

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

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

[math] \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]C_{ij}[/math], соответствующую произведению [math]AB[/math].

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

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

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

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

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

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

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

Число блоков определяется как округление снизу корня квадратного из числа доступных процессов.

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

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

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

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

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

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

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

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

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

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

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

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

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

Возьмем квадратные матрицы размерностью [math]n \times n[/math], где [math]500\leqslant n \leqslant 2000[/math], а число процессов будем изменять от [math]2[/math] до [math]128[/math]. В качестве примера приведем время выполнения программы в зависимости от числа процессов при фиксированном размере матриц [math]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]64[/math] процесса.

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

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

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

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

Для компиляции был использован компилятор языка C++ GNU 4.4.6. Использованная реализация MPI: Intel MPI 4.0.3. Опции компилятора: -lm.

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

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

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

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

[2] - ИНТУИТ

[3] - OpenNET