Участник:Kisel dv/DNS алгоритм умножения матриц

Автор: Д.В.Кисель

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

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

Алгоритм DNS(Dekel, Nassimi, Sahni - назван по фамилиям создателей) является одним из блочных алгоритмов решения задачи [math]\displaystyle{ C = AB }[/math]  перемножения двух матриц. Для данного алгоритма мы не используем какие-либо свойства матриц [math]\displaystyle{ A }[/math]  и [math]\displaystyle{ B }[/math] . DNS основан на блочном разделении данных, его преимуществом является временная сложность [math]\displaystyle{ O(log(n)) }[/math]  при вычислительной сложности [math]\displaystyle{ O(\frac{n^3}{log(n)}) }[/math] . Достигается это за счет 3d-секционирования(англ. partitioning), в отличие от альтернатив, использующих 2d-секционирование(например, алгоритм Кэннона). Тем самым, алгоритм DNS представим в виде куба, где матрицы [math]\displaystyle{ A }[/math] , [math]\displaystyle{ B }[/math]  и матрица [math]\displaystyle{ C }[/math]  ортогональны друг другу. Далее данное описание будет рассмотрено подробнее.

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

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

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

Используемые формулы:

[math]\displaystyle{ C_{ij} = \sum_{k=0}^{m-1} A_{ik} B_{kj} \quad i \in [0, m-1], \quad j \in [0, m-1] }[/math] , где [math]\displaystyle{ m \times m }[/math]  - количество блоков матрицы.

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

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

[math]\displaystyle{ \sum_{k=0}^{m-1} A_{ik} B_{kj} \quad i \in [0, m-1], \quad i \in [0, m-1] }[/math] 

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

В алгоритме используется:

[math]\displaystyle{ \cdot }[/math]  произведение блоков матриц [math]\displaystyle{ A }[/math]  и [math]\displaystyle{ B }[/math] 

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

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

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

[math]\displaystyle{ C_{ij} = \sum_{k=0}^{m-1} A_{ik} B_{kj} \quad i \in [0, m-1], \quad j \in [0, m-1] }[/math] 

В результате получаем матрицу [math]\displaystyle{ C }[/math] , равную произведению матриц [math]\displaystyle{ A }[/math]  и [math]\displaystyle{ B }[/math] .

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

Рассматриваем матрицы размером [math]\displaystyle{ n \times n }[/math] , разбитые на [math]\displaystyle{ m \times m }[/math]  блоков размера [math]\displaystyle{ \frac{n}{m} }[/math] .

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

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

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

Серым на графе помечен этап, в котором происходит считывание матриц [math]\displaystyle{ A }[/math] и [math]\displaystyle{ B }[/math] и их запись в "смежные грани куба".

Далее как бы проецируем каждый блок матрицы на каждую плоскость куба из процессов(на каждую плоскость, параллельную рассматриваемой матрице, см.рисунок выше), таким образом в каждом из них будет выполняться свое умножение уникальной комбинации блоков("синий" этап).

Зеленым обозначен этап суммирования и вывода результата - матрицы [math]\displaystyle{ C }[/math].

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

[math]\displaystyle{ \cdot }[/math]  Рассмотрим [math]\displaystyle{ m^3 }[/math]  процессов, доступных для умножения двух матриц размером [math]\displaystyle{ n \times n }[/math] , разбитых на [math]\displaystyle{ m \times m }[/math]  блоков. Представим их в виде трехмерного массива [math]\displaystyle{ m \times m \times m }[/math] .

[math]\displaystyle{ \cdot }[/math]  Процессы именованы согласно их расположению в массиве, соответственно, вычисление произведения [math]\displaystyle{ A_{ik}B_{kj} }[/math]  присвоено процессу с координатами [math]\displaystyle{ [i, j, k] (0 \leq i,j,k \lt m) }[/math] .

[math]\displaystyle{ \cdot }[/math]  После того, как каждый процесс завершил умножение, результаты процессов [math]\displaystyle{ [i, j, 0], [i, j, 1], ..., [i, j, m-1] }[/math]  суммируются в [math]\displaystyle{ C_{ij} }[/math] . Суммирование всех [math]\displaystyle{ C_{ij} }[/math]  может выполняться одновременно за [math]\displaystyle{ log(m) }[/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 }[/math], где [math]\displaystyle{ 500\leqslant N \leqslant 2000 }[/math], а число процессов изменялось от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ 128 }[/math]. Для примера посмотрим как изменялось время выполнения программы в зависимости от числа процессов (здесь берутся матрицы размера [math]\displaystyle{ N = 2000 }[/math])

Размер блоков, на которые разбивается матрица, выбирается по оптимальному принципу: целая часть от [math]\displaystyle{ \frac{N}{q^{\frac{1}{3}}} }[/math], где [math]\displaystyle{ N, q }[/math], соответственно, размер матрицы и количество процессов.

Число процессов	Время (с)
128	0.976
64	2.245
32	6.012
16	18.012
8	18.153
4	162.214

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

Проиллюстрируем данные графиком:

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

Размер матриц	Время (с)
2000	2.245
1500	0.814
1000	0.232
750	0.094
500	0.027

Проиллюстрируем данные графиком:

3d-график:

Очевидно, что система хорошо масштабируема.

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

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

Для компиляции был использован компилятор языка C++ GNU 4.4.6, реализация MPI: Intel MPI 4.0.3. При компиляции использована опция -lm (подключение библиотеки math).

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

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

DNS алгоритм не очень распространен в Интернете, сравнительно просто отыскать лишь текстовое описание данного способа умножения матриц. Чуть более трудной, однако выполнимой задачей будет нахождение его реализации, например, на github.