Министерство образования и науки
Государственное образовательное учреждение высшего профессионального образования «Нижегородский государственный университет
им. Н.И. Лобачевского»
Факультет вычислительной математики и кибернетики
Кафедра: Математического обеспечения ЭВМ
Направление: Информационные технологии
Отчет по лабораторной работе
Тема:
«Умножение разреженных матриц»
Выполнили:
студенты группы 85М21
Леденцов Д. К.
Орлов Л. М.
Нижний Новгород
2011
Оглавление
Отчет по лабораторной работе 1
Введение 3
Постановка задачи 4
Описание алгоритма 5
Параллельная реализация 8
Реализация на CUDA 10
Реализация на TBB 12
Реализация на OpenCL 14
Полученные результаты 16
Заключение 18
Литература 19
Введение
По мере того, как растут производительность и быстродействие вычислительных машин, становится возможным обрабатывать все большего размера матрицы для различного типа задач. Например, обработка матриц необходимо при решении систем линейных уравнений. Но, несмотря на стремительное развитие вычислительной техники, по-прежнему, как и несколько десятков лет назад, основными характеристиками остаются: память, трудоемкость и быстродействие. С ростом порядка матричной задачи растет и стоимость ее решения, становясь решающим фактором.
Разреженной матрицей называется матрица с большим числом ненулевых элементов. Для такой матрицы неэффективно хранить и обрабатывать все элементы. Поэтому существуют специальные форматы хранения данных матриц, позволяющие запоминать только ненулевые элементы. Но для каждого такого нового способа хранения необходимо разрабатывать методы для основных операций: умножение матрицы на число, транспонирование, умножение двух матриц.
Постановка задачи
В ходе выполнения лабораторной работы необходимо реализовать умножение разреженных матриц. При этом матрицы должны храниться в разреженном столбцовом формате. Изначально требуется реализовать последовательный вариант алгоритма умножения на языке С++. Далее необходимо разработать параллельную версию алгоритма и реализовать ее с помощью трех технологий: CUDA, TBB, OpenCL. В завершении работы должны быть приведены полученные результаты и заключение.
Описание алгоритма
Для хранения разреженных матриц используется так называемый разреженный столбцовый формат, широко известным как CSS (Compressed Column Storage) или CSC (Compressed Sparse Columns). В таком формате матрица хранится в виде трех массивов:
-
Массив значений Value (построчно, сверху вниз);
-
Массив номеров строк Row;
-
Массив индексов начала столбцов ColIndex.
При таком хранении ColIndex[j] указывает на начало j-го столбца в массивах Value и Row, а элементы этого столбца находятся по индексам от ColIndex[j] до ColIndex[j+1]–1 включительно. При этом довольно просто обрабатываются пустые столбцы, для которых ColIndex[j] = ColIndex[j+1]. Для последнего столбца полагается ColIndex[N+1] = NZ, где N - размер матрицы, а NZ - число ненулевых элементов.
Для хранения матрицы N*N с NZ ненулевыми элементами в итоге нам потребуется 8NZ+4NZ+4(N+1) = 12NZ+4N+4 байт, если для хранения значений элементов мы будем использовать 8 байт, а для индексов - 4 байта. Это конечно много меньше N^2 как если бы мы использовали плотный формат.
Ниже представлен пример хранения для матрицы A (6*6).
Рис.1. Столбцовый формат хранения разреженной матрицы
Формат CSS предоставляет быстрый доступ к столбцам, при этом столбцы рассматриваются по порядку, хотя внутри них элементы могут быть, как упорядочены, так и нет. В первом случае достигается быстрый поиск элементов, хотя приходится “платить” за поддержание такой упорядоченности. Во втором же случае ничего поддерживать не надо, но зато поиск осуществляется лишь перебором. Существует также модификация CSS с четырьмя массивами, где последний массив хранит индексы элементов, идущих в конце столбца. Столбцы тогда могут не быть упорядочены, а их перестановка проводится без перепаковки путем изменения индексов. Данная модификация в лабораторной работе не используется.
Для реализации умножения нужно, прежде всего, реализовать транспонирование матрицы. Пусть имеется разреженная матрица A(N*N) в формате CSS. Чтобы получить AT в формате CSS необходимо:
-
сформировать N «целых» и N «вещественных» векторов;
-
в цикле просмотреть все столбцы исходной матрицы, для каждого столбца - все его элементы;
-
если A[i][j] = v, тогда добавить числа i и v в j-ые «целый» и «вещественный» вектор, тем самым в векторах сформируются столбцы транспонированной матрицы;
-
скопировать данные из векторов в CSS-структуру транспонированной матрицы (Row и Value), попутно формируя массив ColIndex.
Теперь можно непосредственно выполнять умножение. Алгоритм умножение разреженных матриц A и B в форматах CSS состоит в следующем:
-
транспонировать матрицу A, т.е. вычислить AT;
-
инициализировать структуру данных для матрицы C = A*B;
-
последовательно перемножить каждый столбец матрицы AT на каждую из строк матрицы B, записывая в C полученные результаты и формируя ее структуру.
При умножении столбцов возникает задача сопоставления с целью выделения пар ненулевых элементов. В простейшем варианте это решается простым перебором для каждого элемента столбца матрицы AT элементов столбца матрицы B до тех пор, пока не будет найден элемент с таким же значением в массиве Row или не закончится строка. Хотя это излишне, ведь вектора упорядочены! Для решения проблемы достаточно:
-
встать на начало обоих векторов (ks= …, ls= …);
-
сравнить текущие элементы AT.Row[ks]и B.Row[ls];
-
если значения совпадают, просуммировать AT.Value[ks] * B.Value[ls] и увеличить оба индекса, в противном случае – увеличить один из индексов, в зависимости от того, какое значение больше.
Этого можно добиться, написав приблизительно следующее:
ks:=AT.ColIndex[i]; kf:=AT.ColIndex[i + 1] -1
ls:= B.ColIndex[j]; lf:= B.ColIndex[j + 1] -1
while ((ks <= kf) && (ls <= lf))
if AT.Row[ks]
else if AT.Row[ks]>B.Row[ls] then ls:=ls+1;
else sum:=sum + AT.Value[ks]*B.Value[ls];
ks:=ks+1;ls:=ls+1;
endif;
endif
endwile
Параллельная реализация
Предположим, что нам необходимо умножить две разреженные матрицы размера N*N, которые хранятся в столбцовом формате.
C = A * B
Для этого сначала нужно транспонировать матрицу A, а затем перемножить столбцы транспонированной матрицы A на столбцы матрицы B. Большую часть времени работы алгоритма занимает умножение матриц. Поэтому данная часть будет распараллелена. Транспонирование будет выполняться последовательно.
Умножение
После того, как выполнили транспонирование матрицы A, имеются следующие исходные данные:
-
массив valueA – ненулевые элементы матрицы Ат;
-
массив rowA – массив номеров строк элементов матрицы Ат;
-
массив colIndexA – массив начала столбцов матрицы Ат;
-
массив valueB – ненулевые элементы матрицы B;
-
массив rowB – массив номеров строк элементов матрицы B;
-
массив colIndexB – массив начала столбцов матрицы B;
Умножение матриц выполняется в два прохода. Количество потоков равняется размеру матриц (N). На первом проходе каждый поток подсчитывает количество ненулевых элементов в соответствующем столбце матрицы C и записывает это значение в элемент colIndexC[i]. После этого один поток пересчитывает массив colIndexC по следующему алгоритму:
int t;
int sum = 0;
for (int i=0; i<=N; i++)
{
t = colIndexC[i];
colIndexC[i] = sum;
sum += t;
}
В результате будет сформирован массив colIndex для матрицы C. Последний элемент данного массива будет содержать количество ненулевых элементов в матрице С.. Далее будет выделена память под массивы: rowC и valueC.
На втором проходе каждый вычисляет ненулевые элементы матрицы С для своего столбца и записывает их в массив valueC а номер строки в массив rowC. После этого матрица C будет сформирована.
Реализация на CUDA
Для хранения структуры разреженной матрицы был создан класс SparseMatrix. Данный класс содержит следующие поля:
-
int N – размер матрицы;
-
int size – количество ненулевых элементов;
-
float* Value – указатель на массив элементов матрицы;
-
int* Row – указатель на массив номеров строк;
-
int* colIndex – указатель на массив начала столбцов.
Также в данном классе реализованы два метода:
-
virtual SparseMatrix Transpose() - транспонировние матрицы;
-
virtual SparseMatrix Multiplication(SparseMatrix* m) – перемножение двух разреженных матриц.
Для реализации алгоритма перемножения разреженных матриц с помощью технологии CUDA был создан класс SparseMatrix_CUDA, который является наследником класса SparseMatrix,и в нем переопределен метод Multiplication(). Данный метод выделяет память на устройстве, и вызывает функции ядра, выполняющие траснпонирование и умножение матриц. Затем полученные данные считываются с устройства и возвращается результирующая матрица. Функции ядра имеют следующий прототип:
__global__ void FillTMatrix(int n, float* value, int* row, int* colIndex, float* tMatrix)
где n – размер матрицы, value – указатель на область памяти, в которой хранится массив значений транспонируемой матрицы, row – указаетель на массив номеров строк, colIndex – указатель на массив начала столбцов, tMatrix – транспонированная матрица. Данная функция заполнеят транспонированную матрицу.
__global__ void FillColIndex(int n, int* colIndex, float* tMatrix)
где n – размер матрицы, colIndex – указатель на массив начала столбцов транспонированной матрицы, tMatrix – элементы транспонированной матрицы. Данная функция записывает в массив colIndex количество ненулевых элементов в каждом столбце.
__global__ void RecalcColIndex(int n, int* colIndex)
где n – размер матрицы, colIndex – указатель на массив начала столбцов транспонированной матрицы. Данную функцию выполняет один поток, который пересчитывает массив colIndex.
__global__ void CalcRowAndValue(int n, float* value, int* row, int* colIndex, float* tMatrix)
где n – размер матрицы, value – указатель на область памяти, в которой хранится массив значений транспонируемой матрицы, row – указаетель на массив номеров строк, colIndex – указатель на массив начала столбцов, tMatrix – транспонированная матрица. Данная функция заполняет массивы row и value у транспонированной матрицы.
__global__ void MultiplicationGPU(int n, float* valueA, int* rowA, int* colIndexA, float* valueB, int* rowB, int* colIndexB, float* resMatrix),
где n – размер матриц, valueA, rowA, colIndexA – массивы матрицы A, valueB, rowB, colIndexB – масивы матрицы B, resMatrix – полученная матрица.
В данных функциях вычисления выполняются по описанным выше алгоритмам. Для выполнения ядер транспонирования на устройстве создается одномерное пространство индексов. И номер потока определяется следующим образом:
int glIdx = blockIdx.x * blockDim.x + threadIdx.x;
Для выполнения ядра умножения на устройстве создается двумерное пространство индексов и номер потока определяется следующим образом:
int j = blockIdx.x * blockDim.x + threadIdx.x;
int i = blockIdx.y * blockDim.y + threadIdx.y;
Реализация на TBB
Для реализации алгоритма перемножения разреженных матриц с помощью технологии TBB был создан класс SparseMatrix_TBB, который является наследником класса SparseMatrix,и в нем был переопределен метод Multiplication().
Для распараллеливания алгоритма умножения матриц с помощью библиотеки TBB был создан класс Multiplicator. Данный класс содержит следующие данные:
-
vector* rows - массив векторов типа int из библиотеки STL. Данные вектора будут содержать номера строк ненулевых элементов матрицы C.
-
vector* values - массив векторов типа float из библиотеки STL. Данные вектора будут содержать ненулевые элементы матрицы C.
-
int *col_index – указатель на массив colIndex матрицы C.
У класса Multiplicator перегружен оператор:
void operator()(const blocked_range &r) const ;
В данном методе заполняются массивы векторов rows, values, и в массив colIndex записывается количество ненулевых элементов для каждого столба. Реализация данного метода представлена ниже.
void operator()(const blocked_range& r) const
{
int begin = r.begin();
int end = r.end();
int N = A.N;
int i, j, k;
int *temp = new int[N];
for (i = begin; i < end; i++)
{
memset(temp, -1, N * sizeof(int));
int ind1 = A.ColIndex[i], ind2 = A.ColIndex[i + 1];
for (j = ind1; j < ind2; j++)
{
int row = A.Row[j];
temp[row] = j;
}
for (j = 0; j < N; j++)
{
double sum = 0;
int ind3 = B.ColIndex[j], ind4 = B.ColIndex[j + 1];
for (k = ind3; k < ind4; k++)
{
int brow = B.Row[k];
int aind = temp[brow];
if (aind != -1)
sum += A.Value[aind] * B.Value[k];
}
if (fabs(sum) > 0.000001)
{
rows[i].push_back(j);
values[i].push_back(sum);
col_index[i]++;
}
}
}
delete [] temp;
}
Таким образом, для выполнения умножения разреженных матриц с помощью технологии TBB необходимо сделать:
-
Создать массивы rows и values размера N из векторов типа int и float.
-
Вызвать метод parallel_for, в который передать экземпляр класса Multiplicator. После этого массив colIndex результируещей матрицы будет содержать количество ненулевых элементов для каждого столбца.
-
Пересчитать массив colIndex.
-
Создать матрицу и передать в нее данные из массивов: colIndex, rows, values.
Реализация на OpenCL
Для реализации алгоритма перемножения разреженных матриц с помощью технологии OpenCL был создан класс SparseMatrix_OpenCL, который является наследником класса SparseMatrix, и в нем переопределен метод Multiplication(). В данном методе выбирается платформа для вычислений, устройство, создается контекст, очередь команд, ядра, буферы памяти. Сначала выполняется последовательный алгоритм транспонирования матрицы A. Затем выполняются функции ядра для умножения матриц.
Для выполнения умножения матриц необходимы следующие буферы памяти:
-
cl_mem maValue – массив элементов матрицы Ат;
-
cl_mem maRow – массив номеров строк элементов матрицы Ат;
-
cl_mem maColIndex – массив начала столбцов матрицы Ат;
-
cl_mem mbValue – массив элементов матрицы B;
-
cl_mem mbRow - массив номеров строк элементов матрицы B;
-
cl_mem mbColIndex - массив начала столбцов матрицы B;
-
cl_mem mсValue – массив элементов матрицы С;
-
cl_mem mсRow - массив номеров строк элементов матрицы С;
-
cl_mem mсColIndex - массив начала столбцов матрицы С;
Первое ядро __kernel void Mult_CalcColIndex выполняет подсчет количества ненулевых элементов для столбцов матрицы С. Эти данные записываются в буфер mcColIndex. Код ядра приведен ниже:
__kernel void Mult_CalcColIndex(const int n,
__global float* valueA,
__global int* rowA,
__global int* colIndexA,
__global float* valueB,
__global int* rowB,
__global int* colIndexB,
__global int* colIndexC)
{
int j = get_global_id(0) ;
if (j >= n) return;
int ks, ls, kf, lf;
int countElem = 0;
for (int i=0; i
{
ks = colIndexA[i];
kf = colIndexA[i+1];
ls = colIndexB[j];
lf = colIndexB[j+1];
bool isFound = false;
while ((ks < kf) && (ls < lf))
{
if (rowA[ks] < rowB[ls]) ks++;
else
if (rowA[ks] > rowB[ls]) ls++;
else
{
isFound = true;
break;
}
}
if (isFound) countElem++;
}
colIndexC[j] = countElem;
}
Затем массив colIndexC считывается с устройства и пересчитывается на хосте и запускается второе ядро __kernel void Mult_CalcValue. Оно вычисляет массивы valueC и rowC.
__kernel void Mult_CalcValue(const int n,
__global float* valueA,
__global int* rowA,
__global int* colIndexA,
__global float* valueB,
__global int* rowB,
__global int* colIndexB,
__global float* valueC,
__global int* rowC,
__global int* colIndexC)
{
int j = get_global_id(0) ;
if (j >= n) return;
int ks, ls, kf, lf;
int idxCurElem = colIndexC[j];
for (int i=0; i
{
ks = colIndexA[i];
kf = colIndexA[i+1];
ls = colIndexB[j];
lf = colIndexB[j+1];
bool isFound = false;
float sum = 0;
while ((ks < kf) && (ls < lf))
{
if (rowA[ks] < rowB[ls]) ks++;
else
if (rowA[ks] > rowB[ls]) ls++;
else
{
isFound = true;
sum += valueA[ks] * valueB[ls];
ks++;
}
}
if (isFound)
{
valueC[idxCurElem] = sum;
rowC[idxCurElem] = i;
idxCurElem++;
}
}
}
После этого полученные данные считываются с устройства, и создается матрица C.
Полученные результаты
Для сравнения последовательной и параллельных реализаций алгоритма умножения разреженных матриц была проведена серия экспериментов. Тестирование алгоритма CUDA-версии проводилось на системе со следующими характеристиками:
-
CPU – Intel(R) Core(TM) Quad 2.4 GHz
-
RAM – 4 GB
-
GPU – NVidia GeForce 8600GT
-
OS – Windows Seven
В таблице приведено время работы алгоритмов и полученное ускорение в зависимости от размера матриц. Необходимо отметить, что время работы CUDA-версии включает время, затрачиваемое на выделение памяти и копирование данных.
Размер матриц
|
Последовательная версия (сек.)
|
CUDA-версия (сек.)
|
Ускорение
|
128*128
|
0,0172
|
0,155
|
0,11
|
256*256
|
0,132
|
0,167
|
0,79
|
512*512
|
1,037
|
0,29
|
3,58
|
1024*1024
|
8,28
|
1,4
|
5,91
|
Таблица . Результаты CUDA-версии
Тестирование алгоритмов TBB-версии и OpenCL-версии проводилось на системе со следующими характеристиками:
-
CPU – AMD Turion(tm) 64 X2 2.1 GHz
-
RAM – 3 GB
-
GPU – Radeon HD 3470
-
OS – Windows Vista
Ниже в таблицах приведено время работы алгоритмов и полученное ускорение в зависимости от размера матриц. Ядро алгоритма OpenCL-версии выполнялось на центральном процессоре.
Размер матриц
|
Последовательная версия (сек.)
|
OpenCL-версия
|
TBB-версия
|
128*128
|
0,018
|
0,552
|
0,005
|
256*256
|
0,148
|
0,596
|
0,03
|
512*512
|
1,182
|
1,001
|
0,202
|
1024*1024
|
9,046
|
4,138
|
1,624
|
2048*2048
|
56,017
|
26,283
|
18,108
|
Таблица . Время работы TBB-версии и OpenCL-версии
Размер матриц
|
Ускорение OpenCL-версии
|
Ускорение TBB-версии
|
128*128
|
0,033
|
3,6
|
256*256
|
0,248
|
4,93
|
512*512
|
1,18
|
5,85
|
1024*1024
|
2,19
|
5,57
|
2048*2048
|
2,13
|
3,09
|
Таблица . Ускорение TBB-версии и OpenCL-версии
Заключение
В ходе выполнения лабораторной работы был успешно реализован последовательный вариант алгоритма умножения разреженных матриц на языке С++. Также разработаны параллельные версии данного алгоритма, использующие технологии: CUDA, TBB, OpenCL.
Проведенные эксперименты для CUDA-версии показали, что параллельная версия алгоритма начинает работать быстрее при размере матрицы больше 512. Это связано с тем, что вычисление на видеокарте требует накладных расходов из-за выделения памяти и копирования данных. Таким образом, данный алгоритм эффективен для матриц больших размеров.
Алгоритмы TBB-версии и OpenCL-версии выполнялись на центральном процессоре. Как видно из полученных результатов, при небольшом размере матриц Open-CL версия работает медленнее последовательной версии. Это связано с накладными расходами на создание объектов библиотеки OpenCL, выделение и копирование памяти, компиляцию ядер. Но с ростом размера матриц ускорение стремится к двум. Время работы TBB-версии алгоритма в несколько раз превосходит время работы последовательного алгоритма. Помимо выигрыша связанного с распараллеливанием, ускорение еще достигается благодаря оптимизации алгоритма. Алгоритм TBB-версии выполняет умножение матриц не в два прохода, а в один.
Литература
-
Писсанецки С. Технология разреженных матриц.: Пер. с англ. - М.: Мир, 1988.
-
Тьюарсон Р. Разреженные матрицы.: Пер. с англ. М.: Мир, 1977.
-
Боресков, А.В. Основы работы с технологией CUDA / А.В. Боресков, А.А. Харламов. – М.: ДМК Пресс, 2011. – 232 с
Интернет – ресурсы:
-
Разреженные матрицы [http://en.wikipedia.org/wiki/Sparse_matrix]
-
Технология CUDA [http://developer.nvidia.com/]
-
Технология OpenCL для AMD [http://www.amd.com/us/products/technologies/stream-technology/opencl/Pages/opencl.aspx]
-
Технология TBB [www.threadingbuildingblocks.org]
|
