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

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

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

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

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

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

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

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

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

$\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}$

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

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

$\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}$

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

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

$\begin{align} \sum_{k = 1}^{n} a_{ik} b_{kj} \end{align}$

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

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

$\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}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Объём входных данных: $2n^{2}$

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

Объём выходных данных: $n^{2}$

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

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

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

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

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

Далее рассмотрим зависимость времени выполнения программы от размера матриц. Число процессов также зафиксируем: возьмем $64$ процесса.

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

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

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

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

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

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

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

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

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

[2] - ИНТУИТ

[3] - OpenNET