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

Выполнила: И.В. Близнякова (611 группа).

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

Общее описание алгоритма^[1]

Разрежённая матрица — это матрица с преимущественно нулевыми элементами. В противном случае, если бо́льшая часть элементов матрицы ненулевые, матрица считается плотной.

Среди специалистов нет единства в определении того, какое именно количество ненулевых элементов делает матрицу разрежённой. Разные авторы предлагают различные варианты. Для матрицы порядка n число ненулевых элементов:

есть [math]\displaystyle{ O(n) }[/math]. Такое определение подходит разве что для теоретического анализа асимптотических свойств матричных алгоритмов;
в каждой строке не превышает 10 в типичном случае;
ограничено [math]\displaystyle{ n^{1+\gamma} }[/math], где [math]\displaystyle{ \gamma \lt 1 }[/math].
таково, что для данного алгоритма и вычислительной системы имеет смысл извлекать выгоду из наличия в ней нулей.

Хранение разреженной матрицы

Формат RR(C)O

Рассмотрим сначала формат RR(C)O. Сокращенное название данного формата происходит от английского словосочетания "Row - wise Representation Complete and Ordered" (строчное представление, полное и упорядоченное). В данном формате вместо одного двумерного массива, используются три одномерных. Значения ненулевых элементов матрицы и соответствующие им столбцовые индексы хранятся в этом формате по строкам в двух массивах [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math]. Массив указателей [math]\displaystyle{ IA }[/math], используется для ссылки на компоненты массивов [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math], с которых начинается описание очередной строки. Последняя компонента массива [math]\displaystyle{ IA }[/math] содержит указатель первой свободной компоненты в массивах [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math], т.е. равна числу ненулевых элементов матрицы, увеличенному на единицу. Здесь уместно привести пример.

Рассмотрим матрицу [math]\displaystyle{ A }[/math]:

[math]\displaystyle{ \begin{pmatrix} 0 & 0 & 0 & 2 & 0 \\ 1 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 6 & 0 & 0 & 0 \\ 0 & 0 & 4 & 0 & 0 \\ \end{pmatrix} }[/math],

тогда ее представление в формате RR(C)O будет иметь вид:

  IA = [ 1 2 4 4 5 6 ]
  JA = [ 4 1 3 2 3 ]
  AN = [ 2 1 3 6 4 ]

Т.е. массив [math]\displaystyle{ AN }[/math] содержит все не нулевые элементы исходной матрицы [math]\displaystyle{ A }[/math], массив [math]\displaystyle{ JA }[/math] номер столбца в котором находится соответствующий элемент из [math]\displaystyle{ AN }[/math] и наконец массив [math]\displaystyle{ IA }[/math] содержит номер с которого начинается описание элементов в массивах [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math]. Таким образом информация об элементах 2-ой строки матрицы хранится в элементах с [math]\displaystyle{ IA[2] = 2 }[/math] по [math]\displaystyle{ IA[3] - 1 = 3 }[/math] включительно массивов [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math]. Можно обратить внимание, что [math]\displaystyle{ IA[3] = IA[4] = 4 }[/math], а это означает, что 3-я строка матрицы [math]\displaystyle{ A }[/math] нулевая.

В общем случае описание [math]\displaystyle{ r }[/math]-й строки матрицы A хранится в компонентах с [math]\displaystyle{ IA[r] }[/math] до [math]\displaystyle{ IA[r + 1] - 1 }[/math] включительно массивов [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math]. Если [math]\displaystyle{ IA[r + 1] = IA[r] }[/math], то это означает, что [math]\displaystyle{ r }[/math]-я строка нулевая. Количество элементов в массиве [math]\displaystyle{ IA }[/math] на единицу больше, чем число строк исходной матрицы, а количество элементов в массивах [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math] равно числу ненулевых элементов исходной матрицы.

Данный способ представления называют полным, поскольку представлена вся матрица [math]\displaystyle{ A }[/math], упорядоченным, поскольку элементы каждой строки матрицы [math]\displaystyle{ A }[/math] хранятся в соответствии с возрастанием столбцовых индексов, и строчным, поскольку информация о матрице [math]\displaystyle{ A }[/math] указывается по строкам.

Массивы [math]\displaystyle{ IA }[/math] и [math]\displaystyle{ JA }[/math] представляют портрет (структуру) матрицы [math]\displaystyle{ A }[/math], задаваемый как множество списков смежности ассоциированного с [math]\displaystyle{ A }[/math] графа. Если алгоритм, реализующий какую-либо операцию над разреженными матрицами, разбит на этапы символической обработки, на котором определяется портрет результирующей матрицы, и численной обработки, на котором определяются значения элементов результирующей матрицы, то массивы [math]\displaystyle{ IA }[/math] и [math]\displaystyle{ JA }[/math] заполняются на первом этапе, а массив [math]\displaystyle{ AN }[/math] - на втором.

Формат RR(C)U

Рассмотрим теперь формат RR(C)U.

Сокращенное название данного формата происходит от английского словосочетания "Row - wise Representation Complete and Unordered" (строчное представление, полное, но неупорядоченное). Формат RR(C)U отличается от RR(C)O тем, что в данном случае соблюдается упорядоченность строк, но внутри каждой строки элементы исходных матриц могут храниться в произвольном порядке. Для матрицы [math]\displaystyle{ A }[/math] нашего примера вполне можно было бы использовать и строчное представление, полное, но неупорядоченное такое:

  IA = [ 1 2 4 4 5 6 ]
  JA = [ 4 3 1 2 3 ]
  AN = [ 2 1 3 6 4 ]

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

Замечания

Несколько замечаний по поводу рассмотренных форматов представления:

Очевидно, что представление матрицы в формате RR(C)O так же является и представлением в формате RR(C)U, но не наоборот.
Из представления матрицы в формате RR(C) нельзя получить информацию о точном количестве столбцов исходной матрицы.
Целесообразно (в вопросе экономии памяти) использовать представления RR(C) в случае, если матрица содержит значительное число нулевых элементов.

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

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

Мы рассмотрим умножение разреженной матрицы общего вида, хранимой в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math] на заполненный вектор-столбец.

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

Исходные данные: разреженная матрица общего вида [math]\displaystyle{ A }[/math] с элементами [math]\displaystyle{ a_{ij} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math] и [math]\displaystyle{ j = 1,...,m }[/math]). Заполненный вектор-столбец [math]\displaystyle{ b }[/math] с элементами [math]\displaystyle{ b_{j} }[/math] ([math]\displaystyle{ j =1,...,m }[/math]).

Вычисляемые данные: заполненный вектор-столбец [math]\displaystyle{ c }[/math] с элементами [math]\displaystyle{ c_{i} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math]).

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

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math],

где [math]\displaystyle{ l_{i} }[/math] - количество ненулевых элементов строки [math]\displaystyle{ i }[/math] матрицы [math]\displaystyle{ A }[/math], [math]\displaystyle{ j(k) }[/math] - индекс [math]\displaystyle{ k }[/math]-го ненулевого элемента матрицы [math]\displaystyle{ A }[/math].

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

Вычислительное ядро последовательной версии умножения разреженной матрицы на вектор можно составить из множественных (всего их [math]\displaystyle{ n }[/math]) вычислений скалярных произведений строк матрицы:

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math].

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

Как записано и в описании ядра алгоритма, основную часть метода составляют множественные (всего [math]\displaystyle{ n }[/math]) вычисления сумм:

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math].

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

Далее предполагаем, что разреженная матрица общего вида [math]\displaystyle{ A }[/math] хранится в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math]. Последовательность исполнения метода следующая:

Выполнять для [math]\displaystyle{ i }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n }[/math]

[math]\displaystyle{ c_{i} = 0 }[/math]
[math]\displaystyle{ IAA = IA[i] }[/math]
[math]\displaystyle{ IAB = IA[i + 1] - 1 }[/math]
[math]\displaystyle{ c_{i} = \sum\limits_{k = IAA}^{IAB} AN[k]b_{JA[k]} }[/math].

После этого (если [math]\displaystyle{ i \le n }[/math]) происходит переход к шагу 1 с бо́льшим [math]\displaystyle{ i }[/math].

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

Для умножения разреженной матрицы общего вида, хранимой в форме RR(C)U, размером [math]\displaystyle{ n \times m }[/math] на заполненный вектор [math]\displaystyle{ m \times 1 }[/math] в последовательном (наиболее быстром) варианте требуется:

[math]\displaystyle{ l }[/math] сложений,
[math]\displaystyle{ l }[/math] умножений.

Умножения и сложения составляют основную часть алгоритма.

При классификации по последовательной сложности, таким образом, алгоритм умножения разреженной матрицы на вектор относится к алгоритмам [math]\displaystyle{ O(l) }[/math].

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

Опишем граф алгоритма^[2]^[3]^[4] как аналитически, так и в виде рисунка.

Граф алгоритма состоит из двух групп вершин, расположенных в целочисленных узлах двух областей одной размерности.

Первая группа вершин расположена в двумерной области, соответствующая ей операция вычисляет функцию [math]\displaystyle{ a \cdot b }[/math]. Естественно введённые координаты области таковы:

[math]\displaystyle{ i }[/math] — меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения;
[math]\displaystyle{ j }[/math] — меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ k }[/math], принимая все целочисленные значения,

где [math]\displaystyle{ k = k(i) = IA[i+1]-IA[i] }[/math].

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

[math]\displaystyle{ a }[/math] - элемент входных данных, а именно [math]\displaystyle{ AN[IA_{i}+j-1] }[/math].
[math]\displaystyle{ b }[/math] - элемент входных данных, а именно [math]\displaystyle{ b[IA_{i}+j-1] }[/math].

Результат срабатывания операции является промежуточным данным алгоритма.

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]\displaystyle{ a + b }[/math]. Естественно введённые координаты области таковы:

[math]\displaystyle{ i }[/math] - меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения;
[math]\displaystyle{ j }[/math] - меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ k-1 }[/math], принимая все целочисленные значения,

где [math]\displaystyle{ k = k(i) = IA[i+1]-IA[i] }[/math].

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

[math]\displaystyle{ a }[/math]:
- при [math]\displaystyle{ j = 1 }[/math] - результат срабатывания операции, соответствующей вершине из первой группы, с координатами [math]\displaystyle{ (i, j) }[/math];
- при [math]\displaystyle{ j \gt 1 }[/math] - результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]\displaystyle{ (i, j - 1) }[/math];
[math]\displaystyle{ b }[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатами [math]\displaystyle{ (i, j + 1) }[/math].

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

при [math]\displaystyle{ j \lt k - 1 }[/math] является промежуточным данным алгоритма;
при [math]\displaystyle{ j = k - 1 }[/math] является выходным данным [math]\displaystyle{ c_{i} }[/math].

Описанный граф можно посмотреть на рис.1. Здесь вершины первой группы обозначены красным цветом и отмечены знаком умножения, вершины второй - зелёным цветом и знаком сложения. Вершины, соответствующие входным данным обозначены белым цветом и выходным - синим.

Рисунок 1. Граф алгоритма умножения разреженной матрицы на вектор

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

Для умножения разреженной матрицы общего вида, хранимой в форме RR(C)U, размерности [math]\displaystyle{ n \times m }[/math] на заполненный вектор [math]\displaystyle{ m \times 1 }[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

Не более чем [math]\displaystyle{ m }[/math] сложений и умножений ([math]\displaystyle{ n }[/math] вычислений в каждом из ярусов)

При классификации по высоте ЯПФ, таким образом, алгоритм умножения разреженной матрицы на вектор относится к алгоритмам со сложностью [math]\displaystyle{ O(m) }[/math]. При классификации по ширине ЯПФ его сложность будет [math]\displaystyle{ O(n) }[/math].

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

Входные данные: разреженная матрица общего вида [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n \times m }[/math], хранимая в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math]. Заполненный вектор-столбец [math]\displaystyle{ b }[/math] с элементами [math]\displaystyle{ b_{j} }[/math] размерности [math]\displaystyle{ m \times 1 }[/math].

Объём входных данных: [math]\displaystyle{ 2l+m+n+1 }[/math], где [math]\displaystyle{ l }[/math] - количество ненулевых элементов в матрице [math]\displaystyle{ A }[/math].

Выходные данные: заполненный вектор-столбец [math]\displaystyle{ c }[/math] с элементами [math]\displaystyle{ c_{i} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math]).

Объём выходных данных: [math]\displaystyle{ n }[/math].

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейной.

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

Пусть [math]\displaystyle{ l }[/math] - количество ненулевых элементов матрицы [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n \times m }[/math]. Пусть [math]\displaystyle{ x }[/math] - объем памяти, используемый для хранения значения элемента матрицы, [math]\displaystyle{ y }[/math] - объём памяти, используемый для хранения номера столбца или строки. В таком случае для хранения матрицы в стандартном представлении нам потребуется объем памяти, равный [math]\displaystyle{ x \cdot n \cdot m }[/math], для хранения в формате RR(C) - [math]\displaystyle{ y(n + 1) + (x + y)l }[/math]. Таким образом хранение в формате RR(C) не является эффективным для матриц, в которых [math]\displaystyle{ l \gt \frac{xnm - y(n+1)}{x+y} }[/math].

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

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

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

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

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

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

Существующие реализации:

SciPy [1],

MatLab [2],

Intel MKL [3].

Литература

↑ С.Писсанецки. Технология разреженных матриц. - М.: Мир, 1988.
↑ Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.
↑ Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.
↑ Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.

[1] С.Писсанецки. Технология разреженных матриц. - М.: Мир, 1988.

[2] Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.

[3] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.

[4] Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.

[1]

[2]

[3]

[4]

Участник:Blizn/Хранение ненулевых элементов разреженной матрицы. Умножение разреженной матрицы на вектор.

Содержание

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

Общее описание алгоритма^[1]

Хранение разреженной матрицы

Формат RR(C)O

Формат RR(C)U

Замечания

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

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

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

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

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

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

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

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

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

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

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

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

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

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

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

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

Литература

Навигация

Участник:Blizn/Хранение ненулевых элементов разреженной матрицы. Умножение разреженной матрицы на вектор.

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

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

Хранение разреженной матрицы

Формат RR(C)O

Формат RR(C)U

Замечания

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

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

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

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

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

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

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

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

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

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

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

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

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

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

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

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

Литература

Навигация

Поиск

Общее описание алгоритма^[1]