Участник:Lexaloris/Умножение разреженной матрицы на вектор

Умножение разреженной матрицы на вектор
Последовательный алгоритм
Последовательная сложность	[math]O(k)[/math]
Объём входных данных	[math]k + m[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(m)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/math]

Авторы страницы: A.Д. Новоселов и П.А. Кочетков

Кочетков Павел ответственен за программную реализацию алгоритма, Новоселов Алексей ответственен за оформление статьи.

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

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

1.1.1 Хранение ненулевых элементов разреженной матрицы

1.1.1.1 Разреженный строчной формат

Одной из наиболее ши­роко используемых схем хранения разреженных матриц является разреженный строчный формат. Эта схема предъявляет минимальные требования к памяти и в то же время оказывается очень удобной для умножения разреженной матрицы на вектор. Например, рассмотрим формат хранения разряженной матрицы [math]A[/math]:

Значения нену­левых элементов матрицы и соответствующие столбцовые индексы хранятся в этой схеме по строкам в двух массивах [math]AN[/math] и [math]JA[/math] соответственно. Используется также массив указателей [math]IA[/math], отмечающих позиции массивов [math]AN[/math] и [math]JA[/math], с которых начинается описание очередной строки. Дополнительная компо­нента в [math]IA[/math] содержит указатель первой свободной позиции в [math]JA[/math] и [math]AN[/math].

Таким образом [math]A[/math] представляется в виде:

  IA = [ 1 4 4 6 ]
  JA = [ 3 4 8 6 8 ]
  AN = [ 1 3 5 7 1 ]

Данный способ представления называют полным, поскольку представлена вся матрица [math]A[/math], и упорядоченным, поскольку эле­менты каждой строки хранятся в соответствии с возрастанием столбцовых индексов. Таким образом, это строчное представление, полное и упорядоченное, или сокращенно (RR (С) О).

1.1.1.2 Неупорядоченное представление

Представления разреженных матриц необязательно должны быть упорядочены в том смысле, что, хотя упорядоченность строк поддерживается, внутри каждой строки элементы могут храниться в произвольном порядке. Для матрицы А нашего примера вполне можно было бы использовать и строчное представление, полное, но неупорядоченное (RR (С) U ).

Неупорядоченное представление [math]A[/math]:

  IA = [ 1 4 4 6 ]
  JA = [ 8 3 4 8 6 ]
  AN = [ 5 1 3 1 7 ]

Неупорядоченные представления могут быть очень удобны. Результаты большинства матричных операций получаются не­ упорядоченными, и упорядочение их стоило бы больших затрат машинного времени. В то же время, за немногими исключениями, алгоритмы для разреженных матриц не требуют, чтобы представления были упорядоченными.

1.1.2 Умножение разреженной матрицы на вектор

Пусть [math]N[/math] — число строк матрицы. Для каждой ее строки [math]I[/math] матрицы мы находим с помощью [math]IA[/math] значения первой [math]IAA[/math] и последней [math]IAB[/math] позиций, занимаемых элементами строки [math]I[/math] в массивах [math]JA[/math] и [math]AN[/math]. Затем, чтобы вычислить скалярное произведение строки [math]I[/math] и вектора [math]B[/math], мы просто просматриваем [math]JA[/math] и [math]AN[/math] на отрезке от [math]IAA[/math] до [math]IAB[/math]: каждое значение, хранимое в [math]JA[/math], есть столбцовый индекс и используется для извлечения из массива [math]B[/math] элемента, который должен быть умножен на соответствующее число из [math]AN[/math]. Результат каж­дого умножения прибавляется к [math]C(I)[/math].

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

Исходные данные:

[math]IA, JA, AN[/math] - заданная матрица в форме (RR (С) U);

[math]B[/math] - заданный заполненный вектор;

[math]N[/math] - число строк матрицы.

Выход: [math]C[/math] вектор-произведение размерности [math]N[/math].

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

[math] \begin{align} & IAA_{i} = IA(i), \quad i \in [1, N], \\ & IAB_{i} = IA(i + 1) - 1, \quad i \in [1, N], \\ & c_{i} = \sum\limits_{j = IAA_{i}}^{IAB_{i}} AN(j)B(JA(j)), \quad i \in [1, N] \\ \end{align} [/math]

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

Вычислительным ядром, т.е. той частью алгоритма, на которую приходится основное время его работы, является вычисление значения [math]i[/math]-го элемента [math]c_{i}[/math] вектора-произведения, т.е произведения строки [math]I[/math] матрицы [math]A[/math] и вектора [math]B[/math] по формуле:

[math] \begin{align} & c_{i} = \sum\limits_{j = IAA_{i}}^{IAB_{i}} AN(j)B(JA(j)), \quad i \in [1, N] \\ \end{align} [/math]

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

Основу алгоритма составляет вычисление значения [math]i[/math]-го элемента [math]c_{i}[/math] вектора-произведения:

[math] \begin{align} & c_{i} = \sum\limits_{j = IAA_{i}}^{IAB_{i}} AN(j)B(JA(j)), \quad i \in [1, N] \\ \end{align} [/math]

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

Последовательность исполнения метода следующая:

Далее для всех [math]i[/math] от [math]1[/math] до [math]N[/math] по нарастанию выполняются:

1. [math] c_{i} = 0; IAA = IA(i); IAB = IA(i + 1 ) - 1 [/math]

После этого, если [math](IAB \le IAA)[/math]:

2. Для всех [math]j[/math] от [math]IAA[/math] до [math]IAB[/math] выполняется:

[math]c_{i} = \sum\limits_{j = IAA_{i}}^{IAB_{i}} AN(j)B(JA(j))[/math]

Псевдокод алгоритма:

   FOR I = 1, N                          (1)
      U = 0.                             (2) 
      IAA = IA(I)                        (3) 
      IAB = IA(I + 1 ) - 1               (4)
      IF NOT(IAB.LT.IAA)                 (5)
          FOR J = IAA, IAB               (6)
          C(I) = U + AN(J)*B(JA(J))      (7)

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

Для всего алгоритма потребуется выполнить [math]O(k)[/math] операций, как в строке стр.7 псевдокода алгоритма (п.1.5), где [math]k[/math] - число ненулевых эле­ментов матрицы.

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

На рисунке 1 изображен информационный граф алгоритма.

Рисунок 1. Информационный граф алгоритма.

По оси [math]X[/math] для отдельной строки [math]I[/math] матрицы отложены в пределах первой [math]IAA[/math] и последней [math]IAB[/math] позиций ненулевых элементов значения [math]AN[/math] и соответствуюие им по столбцовомым индексам, хранимимым в [math]JA[/math], значения [math]B[/math] элементов. Пары этих значений обозначены желтыми квадратиками [math](in)[/math]. Значения in попарно скалярно перемножаются (зеленый кружочек) и суммируются (синий кружочек) в результат значения элемента вектора-произведения [math]c_{i}[/math], т.е [math]out[/math](красный квадратик) с строковым индексом, соответствующем номеру строки матрицы.

По оси [math]Y[/math] отложена схема вычисление значения каждого элемента вектора-произведения для каждой строки матрицы.

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

Алгоритм обладает возможностью выполняться параллельно, что позволяет значительно ускорить вычисления. Вычисление значения каждого элемента вектора-произведения для каждой строки матрицы и элемента вектора [math]B[/math] можно проводить на отдельных процессорах. Максимальная эффективность вычислений достигается при наличии не менее, чем [math]N[/math] процессоров, т.е отдельный процессор на вычисление значения каждого элемента вектора-произведения. На процессоре выполняется [math]k(i)[/math] последовательных операций умножения и сложения. Поэтому число ярусов для [math]i[/math]-ого процессора равно [math]k(i)[/math].

Для параллельного выполнения алгоритма требуется:

• [math]N[/math] ярусов умножений и сложений,
• в каждом из ярусов [math]k \lt \lt m[/math] операций.

При классификации по высоте ЯПФ, алгоритм имеет линейную сложность. При классификации по высоте ЯПФ также линейную.

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

Входные данные:

• [math]IA, JA, AN[/math] - заданная матрица в форме (RR (С) U) c [math]k[/math] ненулевыми элементами и [math]N[/math] строками, [math]m[/math] столбцами;

• [math]B[/math] - заданный заполненный вектор c [math]m[/math] элементами;

Объем входных данных: [math]k + m[/math]

Выходные данные: [math]C[/math] вектор-произведение размерности [math]N[/math].

Объем выходных данных: [math]N[/math]

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

Алгоритм в рамках выбранной версии полностью детерминирован.

Вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – константа.

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

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

Параллельная реализация алгоритма позволяет до некоторой степени сократить время вычислений в несколько раз. На вычислительном комплексе IBM Blue Gene/P было проведено экспериментальное исследование масштабируемости программной реализации умножения разряженной матрицы на вектор. Результаты исследования представлены на рис.2:

Рисунок 2. График зависимости времени выполнения программы умножения разряженной матрицы на вектор от количества ненулевых элементов матрицы для данного числа процессоров.

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

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

Алгоритм релизован в составе библиотек:

• SparseLib++[1]
• cuSPARSE [2]
• uBLAS [3]
• Intel MKL [4]

3 Литература

• [1] С. Писсанецки. Технология разреженных матриц. Изд. Мир, 1988.
• [2] В. В. Воеводин, Вл.В. Воеводин. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002.