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

Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Evgeny Mortikov и ASA.

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

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

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

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

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

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

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

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

1.1.1.1 Формат RR(C)O

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

Рассмотрим матрицу $A$:

$\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}$,

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

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

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

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

Данный способ представления называют полным, поскольку представлена вся матрица $A$, упорядоченным, поскольку элементы каждой строки матрицы $A$ хранятся в соответствии с возрастанием столбцовых индексов, и строчным, поскольку информация о матрице $A$ указывается по строкам.

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

1.1.1.2 Формат RR(C)U

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

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

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

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

1.1.1.3 Замечания

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

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

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

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

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

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

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

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

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

$c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)}$,

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

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

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

$c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)}$.

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

Как записано и в описании ядра алгоритма, основную часть метода составляют множественные (всего $n$) вычисления сумм:

$c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)}$,

которые могут вычисляться в произвольном порядке.

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

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

Выполнять для $i$ от $1$ до $n$

1. $c_{i} = 0$
2. $IAA = IA[i]$
3. $IAB = IA[i + 1] - 1$
4. $c_{i} = \sum\limits_{k = IAA}^{IAB} AN[k]b_{JA[k]}$.

После этого (если $i \le n$) происходит переход к шагу 1 с бо́льшим $i$.

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

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

• $O(l)$ сложений,
• $O(l)$ умножений.

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

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

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

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

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

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

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

где $k = k(i) = IA[i+1]-IA[i]$.

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

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

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

Вторая группа вершин расположена в двумерной области, соответствующая ей операция $a + b$. Естественно введённые координаты области таковы:

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

где $k = k(i) = IA[i+1]-IA[i]$.

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

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

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

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

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

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

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

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

• Не более чем $m$ сложений и умножений ($n$ вычислений в каждом из ярусов)

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

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

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

Объём входных данных: $2l+m+n+1$, где $l$ - количество ненулевых элементов в матрице $A$.

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

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

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

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

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

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

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

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

Параллельная реализация алгоритма заключается в следующем:

1. Каждый процессор считывает входные данные из файла, зависящего от его номера. Данные имеют вид: вектор $b$ и набор определённых строк матрицы $A$ в формате RR(C)U в зависимости от номера процессора и их числа;
2. Процессоры проводят вычисления значений вектора $c$ в зависимости от имеющихся данных;
3. Все процессоры отправляют результаты вычисления процессору, выбранному главным, и главный процессор выводит результат вычислений.

развернуть

Текст использованной в экспериментах реализации (язык С++):

#include <mpi.h>
#include <memory>
#include <iostream>
#include <vector>
#include <sstream>
#include <fstream>
#include <string>
#include <iterator>

typedef std::shared_ptr<std::vector<std::vector<int> > > MyData;

MyData ReadData(const std::string& data_path) {
  std::ifstream file_stream;
  MyData retval(new std::vector<std::vector<int> >());
  file_stream.open(data_path);

  while (!file_stream.eof()) {
    std::string str;
    std::getline(file_stream, str);

    std::istringstream iss(str);
    auto temp_1 = (std::istream_iterator<std::string>(iss));
    auto temp_2 =  (std::istream_iterator<std::string>());

    std::vector<std::string> coordinates{ temp_1, temp_2};
retval->push_back(std::vector<int>());
    for (auto& e : coordinates) {
      retval->back().push_back(std::stoi(e));
    }
  }

  return retval;
}

// input: vector b, vector ia, vector ja, vector an
int main(int argc, char** argv){
	MPI_Init(&argc, &argv);
	int pid, np;
	MPI_Status *s = new MPI_Status;
	MPI_Comm_rank(MPI_COMM_WORLD, &pid);
	MPI_Comm_size(MPI_COMM_WORLD, &np);
	if (argc == 2){
		int n;
		int indx = 0;
		int stp;
                double strttime = MPI_Wtime();
		std::vector<int> b;
		std::vector<int> ia;
		std::vector<int> ja;
		std::vector<int> an;
                auto data = ReadData("data_" + std::to_string(static_cast<long long>(pid)) + ".txt");
	        b = (*data)[0];
		ia = (*data)[1];
		ja = (*data)[2];
		an = (*data)[3];
		n = (*data)[0].size();
		std::vector<int> c; 
                stp = n / np + ((n % np > pid) ? 1 : 0);
		if (pid == 0){
			c.resize(n, 0);
		}
		else{
			c.resize(stp, 0);
		}
		for (int i = 0; i < stp; i++){
			int strt = ia[i] - ia[0];
                        int lstr = ia[i + 1] - ia[i];
			for (int j = 0; j < lstr; j++){
				c[i] += an[strt + j] * b[ja[strt + j]];
			}
		}
		if (pid == 0){
			for (int i = 1; i < np; i++){
				indx += n / np + ((n % np > i - 1) ? 1 : 0);
				stp = n / np + ((n % np > i) ? 1 : 0);
				MPI_Recv(&c[indx], stp, MPI_INT, i, 0, MPI_COMM_WORLD, s);
			}
		}
		else{
			MPI_Send(&c[0], stp, MPI_INT, 0, 0, MPI_COMM_WORLD);
		}
		if (pid == 0){
			for (int e : c) {
                std::cout << e << "\n";
                        }
                        double endtime = MPI_Wtime();
                        std::cout << "Success! Time = " << endtime - strttime << "\n";
		}
	}
	MPI_Finalize();
	return 0;
}

развернуть

Для формирования входных данных использовался следующий алгоритм (язык Python):

from scipy import sparse as sp
from numpy.random import randint
from os import mkdir
from os.path import join as jn

n = 10000
np = [1, 2, 4, 8, 16, 32, 64, 128]

a = sp.csr_matrix(randint(0, 2, (n, n)))
b = randint(1, 5, n)
i = 0

for j in np:
  mkdir('data_' + str(j))
  indx = 0
  for i in range(0, j):
    stp = n / j + (1 if (n % j > i) else 0)
    ael = a.indptr[indx + stp] - a.indptr[indx]
    strt = a.indptr[indx]
    with open(jn('data_' + str(j), 'data_' + str(i) + '.txt'), 'w') as fout:
    
      for e in b:
        fout.write('{0} '.format(e))
      fout.write('\n')
	
      s = a.indptr[indx:indx + stp + 1]
      for e in s:
        fout.write('{0} '.format(e))
      fout.write('\n')
	
      s = a.indices[strt:strt + ael]
      for e in s:
        fout.write('{0} '.format(e))
      fout.write('\n')
    
      s = a.data[strt:strt + ael]
      for e in s:
        fout.write('{0} '.format(e))
      fout.write('\n')
  
    indx = indx + n / j + (1 if (n % j > i) else 0)

Проведём исследование масштабируемости параллельной реализации алгоритма умножения разреженной матрицы на вектор. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Входные данные формировались на персональном компьютере в силу того, что интерпретатор на суперкомпьютере "Ломоносов" не удовлетворял нужным требованиям, и потом отправлялись на суперкомпьютер для последующей обработки.

Программа была собрана на узле компиляции суперкомпьютера "Ломоносов", с помощью компилятора intel/15.0.090 и компилятора openmpi/1.6.5-icc, на языке С++11.

развернуть

Сборка производилась с помощью CMakeLists.txt, его содержание прилагается:

cmake_minimum_required(VERSION 2.8)

project(task1)

set(CMAKE_CXX_STANDARD 11)

set(CXX_STANDARD_REQUIRED)

set(CMAKE_C_COMPILER mpicc)
set(CMAKE_CXX_COMPILER mpicxx)

find_package(MPI REQUIRED)

include_directories(${MPI_INCLUDE_PATH})

SET(CMAKE_CXX_COMPILER mpicc)

if(MPI_COMPILE_FLAGS)
  set_target_properties(task1 PROPERTIES
    COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
endif()

set(CMAKE_CXX_FLAGS "-std=c++11 -fpermissive " ${MPI_COMPILE_FLAGS})

set(SOURCE_EXE task1.cpp)


add_executable(task1 ${SOURCE_EXE})

target_link_libraries(task1)

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

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров [1 : 128] по степеням 2;
• размер матрицы [1000 : 10000] с шагом 1000.

Распределение ненулевых и нулевых коэффициентов: 1:1. Однако стоит обратить внимание, что матрица генерировалась случайным образом, а значит, это лишь распределение, а не точное значение.

Рисунок 2. Время работы программы в зависимости от количества процессоров и размера матрицы

Можно заметить, что при небольшом количестве входных данных, и большом количестве процессоров эффективность распараллеливания падает до нуля, поскольку время вычисления элементов вектора $c$ становится несоизмеримо малым по сравнению с временем, затраченным на пересылку в пункте (3).

Рисунок 4. Эффективность распараллеливания в зависимости от количества процессоров и размера матрицы

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

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

SciPy [1],

MatLab [2],

Intel MKL [3].