Algorithm level

Difference between revisions of "Dense matrix multiplication (serial version for real matrices)"

From Algowiki
Jump to navigation Jump to search
[unchecked revision][checked revision]
 
(12 intermediate revisions by 2 users not shown)
Line 1: Line 1:
Основные авторы описания: [[Участник:Frolov|А.В.Фролов]], [[Участник:VadimVV|Вад.В.Воеводин]] ([[#Описание локальности данных и вычислений|раздел 2.2]]), [[Участник:Teplov|А.М.Теплов]] (раздел [[#Масштабируемость алгоритма и его реализации|2.4]])
+
{{level-a}}
  
== Описание свойств и структуры алгоритма ==
+
Primary authors of this description: [[:ru:Участник:Frolov|A.V.Frolov]].
  
=== Общее описание алгоритма ===
+
== Properties and structure of the algorithm ==
  
'''Перемножение матриц''' - одна из базовых задач в алгоритмах линейной алгебры, широко применяется в большом количестве разных методов.
+
=== General description of the algorithm ===
Здесь мы рассмотрим умножение <math>С = AВ</math>&nbsp; плотных неособенных матриц (последовательный вещественный вариант), то есть тот вариант, где никак не используются ни специальный вид матрицы, ни ассоциативные свойства операции сложения.
 
  
=== Математическое описание ===
+
'''Matrix multiplication''' is one of the basic procedures in algorithmic linear algebra, which is widely used in a number of different methods. Here, we consider the product <math>C = AB</math>&nbsp; of dense real matrices (serial real version), that is, the version in which neither a special form of matrices nor the associative properties of the addition operation are used<ref>Voevodin V.V., Kuznetsov Yu.A. Matrices and computations, Moscow: Nauka, 1984.</ref>.
  
Исходные данные: плотная матрица <math>A</math> (элементы <math>a_{ij}</math>), плотная матрица <math>B</math> (элементы <math>b_{ij}</math>).
+
=== Mathematical description of the algorithm ===
  
Вычисляемые данные: плотная матрица <math>C</math> (элементы <math>c_{ij}</math>).
+
Input data: dense matrix <math>A</math> of size <math>m</math>-by-<math>n</math> (with entries <math>a_{ij}</math>), dense matrix <math>B</math> of size <math>n</math>-by-<math>l</math> (with entries <math>b_{ij}</math>).
  
Формулы метода:
+
Output data: dense matrix <math>C</math> (with entries <math>c_{ij}</math>).
 +
 
 +
Formulas of the method:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
Line 21: Line 22:
 
</math>
 
</math>
  
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
+
There exists also a block version of the method; however, this description treats only the pointwise version.  
  
=== Вычислительное ядро алгоритма ===
+
=== Computational kernel of the algorithm ===
  
Вычислительное ядро перемножения плотных неособенных матриц можно составить из множественных (всего их <math>l</math>) вычислений умножения матрицы <math>A</math> на столбцы матрицы <math>B</math>, или (при более детальном рассмотрении), из множественных (всего их <math>ml</math>) скалярных произведений строк матрицы <math>A</math> на столбцы матрицы <math>B</math>:
+
The computational kernel of the dense matrix multiplication can be compiled of repeated products of <math>A</math> with the columns of <math>B</math> (there are on the whole <math>l</math> such products) or (upon more detailed inspection) of repeated dot products of the rows of <math>A</math> and the columns of <math>B</math> (there are on the whole <math>ml</math> such products):
  
 
:<math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math>  
 
:<math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math>  
  
в режиме накопления или без него, в зависимости от требований задачи.
+
Depending on the problem requirements, these sums are calculated with or without using the accumulation mode.  
  
=== Макроструктура алгоритма ===
+
=== Macro structure of the algorithm ===
  
Как уже записано в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть умножения матрицы на вектор составляют множественные (всего <math>ml</math>) вычисления скалярных произведений строк матрицы <math>A</math> на столбцы матрицы <math>B</math>
+
As already noted in [[#Computational kernel of the algorithm|computational kernel of the algorithm description]], the basic part of the dense matrix multiplication is compiled of repeated dot products of the rows of <math>A</math> and the columns of <math>B</math> (there are on the whole <math>ml</math> such products):
  
 
:<math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math>  
 
:<math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math>  
  
в режиме накопления или без него.
+
These calculations are performed with or without using the accumulation mode.  
  
=== Описание схемы реализации последовательного алгоритма ===
+
=== Implementation scheme of the serial algorithm ===
  
Для всех <math>i</math> от <math>1</math> до <math>m</math>  и для всех <math>j</math> от <math>1</math> до <math>l</math> выполняются
+
For all <math>i</math> from <math>1</math> to <math>m</math>  and for all <math>j</math> from <math>1</math> to <math>l</math>, do
  
 
:<math>c_{ij} = \sum_{k = 1}^{n} a_{ik} b_{kj}</math>
 
:<math>c_{ij} = \sum_{k = 1}^{n} a_{ik} b_{kj}</math>
  
Особо отметим, что вычисления сумм вида <math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math> производят в режиме накопления прибавлением к текущему (временному) значению вычисляемого элемента матрицы <math>c_{ij}</math> произведений <math>a_{ik} b_{kj}</math> для <math>k</math> от <math>1</math> до <math>n</math>, '''c возрастанием''' <math>k</math>, вначале все элементы инициализируются нулями. При суммировании "по убыванию" общая схема принципиально не отличается и потому нами не рассматривается. Другие порядки выполнения суммирования приводят к изменению параллельных свойств алгоритма и будут рассматриваться нами в отдельных описаниях.
+
We emphasize that sums of the form <math>\sum_{k = 1}^{n} a_{ik} b_{kj}</math> are calculated in the accumulation mode by adding the products <math>a_{ik} b_{kj}</math> to the current (temporary) value of <math>c_{ij}</math>.The index <math>k</math> '''increases''' from <math>1</math> to <math>n</math>. All the sums are initialized to zero. The general scheme is virtually the same if summations are performed in the decreasing order of indices; therefore, this case is not considered. Other orders of summation change the parallel properties of the algorithm and are considered in separate descriptions.
 +
 
 +
=== Serial complexity of the algorithm ===
 +
 
 +
The multiplication of two square matrices of order <math>n</math> (that is, <math>m=n=l</math>) in the (fastest) serial version requires
 +
 
 +
* <math>n^3</math> multiplications and the same number of additions.
 +
 
 +
The multiplication of an <math>m</math>-by-<math>n</math> matrix by an <math>n</math>-by-<math>l</math> matrix in the (fastest) serial version requires
 +
 
 +
* <math>mnl</math> multiplications and the same number of additions.
 +
 
 +
The use of accumulation requires that multiplications and additions be done in the double precision mode (or such procedures as the Fortran function DPROD be used). This increases the computation time for performing the matrix multiplication.
 +
 
 +
In terms of serial complexity, the dense matrix multiplication is qualified as a ''cubic complexity'' algorithm (or a ''trilinear complexity'' algorithm if the matrices are rectangular).
 +
 
 +
=== Information graph ===
 +
 
 +
We describe the [[Glossary#Algorithm graph|algorithm graph]] both analytically and graphically.
 +
 
 +
The algorithm graph of multiplying dense matrices consists of a single group of vertices placed at integer nodes of a three-dimensional domain. The corresponding operation is
 +
<math>a+bc</math>.  
  
=== Последовательная сложность алгоритма ===
+
The natural coordinates of this domain are as follows: 
 +
* <math>i</math> varies from <math>1</math> to <math>m</math>, taking all the integer values in this range;
 +
* <math>j</math> varies from <math>1</math> to <math>l</math>, taking all the integer values in this range;
 +
* <math>k</math> varies from <math>1</math> to <math>n</math>, taking all the integer values in this range.
  
Для умножения двух квадратных матриц порядка <math>n</math> (т.е. при <math>m=n=l</math>) в последовательном (наиболее быстром) варианте требуется:
+
The arguments of the operation are as follows:
+
*<math>a</math> is:
* по <math>n^3</math> умножений и сложений.
+
** the constant <math>0</math> if <math>k = 1</math>;
 +
** the result of performing the operation corresponding to the vertex with coordinates <math>i, j, k-1</math> if <math>k > 1</math>;
 +
*<math>b</math> is the element <math>a_{ik}</math> of the ''input data'';
 +
*<math>c</math> is the element <math>b_{kj}</math> of the ''input data'';
  
Для умножения матрицы размером <math>m</math> строк на <math>n</math> столбцов на матрицу размером <math>m</math> строк на <math>n</math> столбцов в последовательном (наиболее быстром) варианте требуется:
+
The result of performing the operation is:
 +
* an intermediate data item if <math>k < n</math>;
 +
* an output data item <math>c_{ij}</math> if <math>k = n</math>.
  
* по <math>mnl</math> умножений и сложений.
+
[[file:Dense mtrx product.png|thumb|center|800px|Figure 1. Multiplication of dense matrices with demonstration of the output data]]
 
При этом использование режима накопления требует совершения умножений и сложений в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения умножения матрицы на вектор.
 
  
При классификации по последовательной сложности, таким образом, алгоритм умножения матриц относится к алгоритмам ''с кубической сложностью'' (в случае неквадратных матриц - с ''трилинейной'').
+
=== Parallelization resource of the algorithm ===
  
=== Информационный граф ===
+
The parallel version of the algorithm for multiplying square matrices of order n requires that the following layers be successively performed:
  
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.
+
* <math>n</math> layers of multiplications and the same numbers of layers for addition (there are <math>n^2</math> operations in each layer).
  
Граф алгоритма умножения плотных матриц состоит из одной группы вершин, расположенной в целочисленных узлах трёхмерной области, соответствующая ей операция  <math>a+bc</math>.
+
The multiplication of an <math>m</math>-by-<math>n</math> matrix by an <math>n</math>-by-<math>l</math> matrix in the (fastest) serial version requires
  
Естественно введённые координаты области таковы:
+
* <math>n</math> layers of multiplications and the same numbers of layers for addition (there are <math>ml</math> operations in each layer).
* <math>i</math> — меняется в диапазоне от <math>1</math> до <math>m</math>, принимая все целочисленные значения;
 
* <math>j</math> — меняется в диапазоне от <math>1</math> до <math>l</math>, принимая все целочисленные значения;
 
* <math>k</math> — меняется в диапазоне от <math>1</math> до <math>n</math>, принимая все целочисленные значения.
 
  
Аргументы операции следующие:
+
The use of accumulation requires that multiplications and subtractions be done in the double precision mode. For the parallel version, this implies that virtually all the intermediate calculations in the algorithm must be performed in double precision if the accumulation option is used. Unlike the situation with the serial version, this results in a certain increase in the required memory size.
*<math>a</math>:
 
** при <math>k = 1</math> константа <math>0</math>;
 
** при <math>k > 1</math> — результат срабатывания операции, соответствующей вершине с координатами <math>i, j, k-1</math>;
 
*<math>b</math> — элемент ''входных данных'', а именно  <math>a_{ik}</math>;
 
*<math>c</math> - элемент ''входных данных'' <math>b_{kj}</math>;
 
  
Результат срабатывания операции является:
+
In terms of the parallel form height, the dense matrix multiplication is qualified as a quadratic complexity algorithm. In terms of the parallel form width, its complexity is also ''quadratic'' (for square matrices) or ''bilinear'' (for general rectangular matrices).
** при <math>k < n</math> - ''промежуточным данным'' алгоритма;
 
** при <math>k = n</math> - выходным данным <math>c_{ij}</math>.
 
  
[[file:Dense mtrx product.png|thumb|center|800px|Умножение плотных матриц с отображением выходных данных]]
+
=== Input and output data of the algorithm ===
  
=== Описание ресурса параллелизма алгоритма ===
+
'''Input data''': matrix <math>A</math> (with entries <math>a_{ij}</math>), matrix <math>B</math> (with entries <math>b_{ij}</math>)).
  
Для алгоритма умножения квадратных матриц порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
+
'''Size of the input data ''': <math>mn+nl</math>
  
* по <math>n</math> ярусов умножений и сложений (в каждом из ярусов — <math>n^2</math> операций).
+
'''Output data ''': matrix <math>C</math> (with entries <math>c_{ij}</math>).
  
Для умножения матрицы размером <math>m</math> строк на <math>n</math> столбцов на матрицу размером <math>n</math> строк на <math>l</math> столбцов в последовательном (наиболее быстром) варианте требуется:
+
'''Size of the output data''': <math>ml</math>
 
* по <math>n</math> ярусов умножений и сложений (в каждом из ярусов — <math>ml</math> операций).
 
  
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
+
=== Properties of the algorithm ===
  
При классификации по высоте ЯПФ, таким образом, алгоритм умножения матрицы на вектор относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''квадратичной'' (для квадратных матриц) или ''билинейной'' (для матриц общего вида).
+
It is clear that, in the case of unlimited resources, the ratio of the serial to parallel complexity is ''quadratic'' or ''bilinear'' (since this is the ratio of the cubic or trilinear complexity to linear one).  
  
=== Описание входных и выходных данных ===
+
The computational power of the dense matrix multiplication, understood as the ratio of the number of operations to the total size of the input and output data, is ''linear''.
  
'''Входные данные''': матрица <math>A</math> (элементы <math>a_{ij}</math>), матрица <math>B</math> (элементы <math>b_{ij}</math>)).
+
The algorithm of dense matrix multiplication is completely determined. We do not consider any other order of performing associative operations in this version.  
  
'''Объём входных данных''': <math>mn+nl</math>
+
== Software implementation of the algorithm ==
  
'''Выходные данные''': матрица <math>C</math> (элементы <math>c_{ij}</math>).
+
=== Implementation peculiarities of the serial algorithm ===
  
'''Объём выходных данных''': <math>ml</math>
+
In its simplest form, the matrix multiplication algorithm in Fortran can be written as follows:
  
=== Свойства алгоритма ===
+
<source lang="fortran">
 +
       
 +
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
 +
</source>
 +
In this case the <math>S</math> variable must be double precision to implement the accumulation mode.
  
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''квадратичным'' или ''билинейным'' (отношение кубической или трилинейной к линейной).
+
=== Possible methods and considerations for parallel implementation of the algorithm ===
 +
=== Run results ===
 +
=== Conclusions for different classes of computer architecture ===
  
При этом вычислительная мощность алгоритма умножения матриц, как отношение числа операций к суммарному объему входных и выходных данных – ''линейно''.
+
== References ==
  
При этом алгоритм умножения матрицы на вектор полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается.
+
<references />
  
 
[[Ru:Перемножение плотных неособенных матриц (последовательный вещественный вариант)]]
 
[[Ru:Перемножение плотных неособенных матриц (последовательный вещественный вариант)]]
  
[[Category:Started articles]]
+
[[Category:Finished articles]]

Latest revision as of 10:55, 8 July 2022


Primary authors of this description: A.V.Frolov.

1 Properties and structure of the algorithm

1.1 General description of the algorithm

Matrix multiplication is one of the basic procedures in algorithmic linear algebra, which is widely used in a number of different methods. Here, we consider the product [math]C = AB[/math]  of dense real matrices (serial real version), that is, the version in which neither a special form of matrices nor the associative properties of the addition operation are used[1].

1.2 Mathematical description of the algorithm

Input data: dense matrix [math]A[/math] of size [math]m[/math]-by-[math]n[/math] (with entries [math]a_{ij}[/math]), dense matrix [math]B[/math] of size [math]n[/math]-by-[math]l[/math] (with entries [math]b_{ij}[/math]).

Output data: dense matrix [math]C[/math] (with entries [math]c_{ij}[/math]).

Formulas of the method:

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

There exists also a block version of the method; however, this description treats only the pointwise version.

1.3 Computational kernel of the algorithm

The computational kernel of the dense matrix multiplication can be compiled of repeated products of [math]A[/math] with the columns of [math]B[/math] (there are on the whole [math]l[/math] such products) or (upon more detailed inspection) of repeated dot products of the rows of [math]A[/math] and the columns of [math]B[/math] (there are on the whole [math]ml[/math] such products):

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

Depending on the problem requirements, these sums are calculated with or without using the accumulation mode.

1.4 Macro structure of the algorithm

As already noted in computational kernel of the algorithm description, the basic part of the dense matrix multiplication is compiled of repeated dot products of the rows of [math]A[/math] and the columns of [math]B[/math] (there are on the whole [math]ml[/math] such products):

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

These calculations are performed with or without using the accumulation mode.

1.5 Implementation scheme of the serial algorithm

For all [math]i[/math] from [math]1[/math] to [math]m[/math] and for all [math]j[/math] from [math]1[/math] to [math]l[/math], do

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

We emphasize that sums of the form [math]\sum_{k = 1}^{n} a_{ik} b_{kj}[/math] are calculated in the accumulation mode by adding the products [math]a_{ik} b_{kj}[/math] to the current (temporary) value of [math]c_{ij}[/math].The index [math]k[/math] increases from [math]1[/math] to [math]n[/math]. All the sums are initialized to zero. The general scheme is virtually the same if summations are performed in the decreasing order of indices; therefore, this case is not considered. Other orders of summation change the parallel properties of the algorithm and are considered in separate descriptions.

1.6 Serial complexity of the algorithm

The multiplication of two square matrices of order [math]n[/math] (that is, [math]m=n=l[/math]) in the (fastest) serial version requires

  • [math]n^3[/math] multiplications and the same number of additions.

The multiplication of an [math]m[/math]-by-[math]n[/math] matrix by an [math]n[/math]-by-[math]l[/math] matrix in the (fastest) serial version requires

  • [math]mnl[/math] multiplications and the same number of additions.

The use of accumulation requires that multiplications and additions be done in the double precision mode (or such procedures as the Fortran function DPROD be used). This increases the computation time for performing the matrix multiplication.

In terms of serial complexity, the dense matrix multiplication is qualified as a cubic complexity algorithm (or a trilinear complexity algorithm if the matrices are rectangular).

1.7 Information graph

We describe the algorithm graph both analytically and graphically.

The algorithm graph of multiplying dense matrices consists of a single group of vertices placed at integer nodes of a three-dimensional domain. The corresponding operation is [math]a+bc[/math].

The natural coordinates of this domain are as follows:

  • [math]i[/math] varies from [math]1[/math] to [math]m[/math], taking all the integer values in this range;
  • [math]j[/math] varies from [math]1[/math] to [math]l[/math], taking all the integer values in this range;
  • [math]k[/math] varies from [math]1[/math] to [math]n[/math], taking all the integer values in this range.

The arguments of the operation are as follows:

  • [math]a[/math] is:
    • the constant [math]0[/math] if [math]k = 1[/math];
    • the result of performing the operation corresponding to the vertex with coordinates [math]i, j, k-1[/math] if [math]k \gt 1[/math];
  • [math]b[/math] is the element [math]a_{ik}[/math] of the input data;
  • [math]c[/math] is the element [math]b_{kj}[/math] of the input data;

The result of performing the operation is:

  • an intermediate data item if [math]k \lt n[/math];
  • an output data item [math]c_{ij}[/math] if [math]k = n[/math].
Figure 1. Multiplication of dense matrices with demonstration of the output data

1.8 Parallelization resource of the algorithm

The parallel version of the algorithm for multiplying square matrices of order n requires that the following layers be successively performed:

  • [math]n[/math] layers of multiplications and the same numbers of layers for addition (there are [math]n^2[/math] operations in each layer).

The multiplication of an [math]m[/math]-by-[math]n[/math] matrix by an [math]n[/math]-by-[math]l[/math] matrix in the (fastest) serial version requires

  • [math]n[/math] layers of multiplications and the same numbers of layers for addition (there are [math]ml[/math] operations in each layer).

The use of accumulation requires that multiplications and subtractions be done in the double precision mode. For the parallel version, this implies that virtually all the intermediate calculations in the algorithm must be performed in double precision if the accumulation option is used. Unlike the situation with the serial version, this results in a certain increase in the required memory size.

In terms of the parallel form height, the dense matrix multiplication is qualified as a quadratic complexity algorithm. In terms of the parallel form width, its complexity is also quadratic (for square matrices) or bilinear (for general rectangular matrices).

1.9 Input and output data of the algorithm

Input data: matrix [math]A[/math] (with entries [math]a_{ij}[/math]), matrix [math]B[/math] (with entries [math]b_{ij}[/math])).

Size of the input data : [math]mn+nl[/math]

Output data : matrix [math]C[/math] (with entries [math]c_{ij}[/math]).

Size of the output data: [math]ml[/math]

1.10 Properties of the algorithm

It is clear that, in the case of unlimited resources, the ratio of the serial to parallel complexity is quadratic or bilinear (since this is the ratio of the cubic or trilinear complexity to linear one).

The computational power of the dense matrix multiplication, understood as the ratio of the number of operations to the total size of the input and output data, is linear.

The algorithm of dense matrix multiplication is completely determined. We do not consider any other order of performing associative operations in this version.

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

In its simplest form, the matrix multiplication algorithm in Fortran can be written as follows:

         
	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

In this case the [math]S[/math] variable must be double precision to implement the accumulation mode.

2.2 Possible methods and considerations for parallel implementation of the algorithm

2.3 Run results

2.4 Conclusions for different classes of computer architecture

3 References

  1. Voevodin V.V., Kuznetsov Yu.A. Matrices and computations, Moscow: Nauka, 1984.