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

Перемножение плотных неособенных матриц
Последовательный алгоритм
Последовательная сложность	$2mnl$
Объём входных данных	$mn+nl$
Объём выходных данных	$ml$
Параллельный алгоритм
Высота ярусно-параллельной формы	$n$
Ширина ярусно-параллельной формы	$2ml$

Основные авторы описания: А.В.Фролов.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$c_{ij} = \sum_{k = 1}^{n} a_{ik} b_{kj}$

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

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

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

• по $n^3$ умножений и сложений.

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

• по $mnl$ умножений и сложений.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Объём входных данных: $mn+nl$

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

Объём выходных данных: $ml$

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

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