Постановка задачи
Матрица — математический объект, записываемый в виде
прямоугольной таблицы элементов кольца или поля, которая представляет собой
совокупность строк и столбцов, на пересечении которых находятся её элементы.
Количество строк и столбцов матрицы задают размер матрицы. Операция умножения
матриц является одной из основных задач матричных вычислений. Умножение матрицы
A размера m×n и матрицы B размера n×k приводит к получению матрицы
С размера m×k, каждый элемент которой определяется скалярным
произведением соответствующей строки матрицы А и столбца матрицы В.
В данной работе рассматриваются только квадратные матрицы,
т. е. m=n=k.
Для этих матриц требуется:
1. Реализовать последовательный алгоритм перемножения
матриц.
2. Реализовать 2 алгоритма ленточного умножения матриц.
3. Провести набор тестов. Сравнить и сопоставить
производительность трех алгоритмов.
Последовательный алгоритм
Последовательный алгоритм умножения матриц представляется
тремя вложенными циклами (матрицу B предварительно
необходимо транспонировать):
void LinearMultiplication(double* A, double* B,
double *C, int
Size)
{
for (int i=0; i < Size; i++)
{
for (int j=0; j < Size;j++)
{
C[i*Size+j]=0;
for
(int k=0; k < Size; k++)
{
C[i*Size+j] +=
A[i*Size+k]*B[j*Size+k];
}
}
}
}
Этот алгоритм является итеративным и ориентирован на
последовательное вычисление строк матрицы С. Предполагается выполнение n·n·n
операций умножения и столько же операций сложения элементов исходных матриц. Количество
выполненных операций имеет порядок O(n3).
Поскольку каждый элемент результирующей матрицы есть
скалярное произведение строки и столбца исходных матриц, то для вычисления всех
элементов матрицы С размером n*n
необходимо выполнить n2·(2n-1)
скалярных операций и затратить время T1 = n2 ·(2n-1)·t, где t есть время выполнения одной
элементарной скалярной операции.
Ленточный алгоритм
Данный алгоритм представляет собой итерационную процедуру.
Исходные матрицы разбиваются на блоки. Количество этих блоков зависит от
количества потоков в программе. На каждой итерации алгоритма все подзадачи
содержат по одной или несколько строк матрицы А и по одному или нескольким
столбцам матрицы В. При выполнении итерации проводится скалярное умножение
содержащихся в подзадачах блоков друг на друга, что приводит к получению
соответствующего блока элементов в вычисляемой матрице С. По завершении
вычислений в конце каждой итерации столбцы матрицы В должны быть переданы между
подзадачами с тем, чтобы в каждой подзадаче оказались новые блоки этой матрицы
для продолжения вычислений матрицы C. На последнем шаге передаются оставшиеся
строки матриц A и B. Таким образом, происходит нахождение матрицы C. Количество
подзадач зависит от количества строк в блоках передаваемых потокам. Например,
если в каждом блоке содержится только по одной строке исходных матриц A и B, то
выделенные базовые подзадачи характеризуются одинаковой вычислительной
трудоемкостью равной n2 (по числу элементов матрицы С).
Когда размер матриц n оказывается больше, чем число
процессоров p, базовые подзадачи можно укрупнить, объединив в рамках одной
подзадачи несколько соседних строк и столбцов перемножаемых матриц. В этом
случае исходная матрица A разбивается на ряд горизонтальных полос, а матрица B
представляется в виде набора вертикальных полос. Размер полос при этом следует
выбрать равным k = n/p (в предположении, что n кратно p), что позволит
по-прежнему обеспечить равномерность распределения вычислительной нагрузки по
процессорам, составляющим многопроцессорную вычислительную систему.
void MatrixMultiplicationMPI(double *&A, double
*&B, double *&C, int &size)
{
int dim =
size;
int i, j, k,
p, ind;
double temp;
MPI_Status Status;
int
ProcPartsize = dim/ProcNum;
int
ProcPartElem = ProcPartsize*dim;
double* bufA =
new double[ProcPartElem];
double* bufB =
new double[ProcPartElem];
double* bufC =
new double[ProcPartElem];
if (ProcRank
== 0)
Transposition(B,size);
MPI_Scatter(A, ProcPartElem,
MPI_DOUBLE, bufA, ProcPartElem, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Scatter(B, ProcPartElem, MPI_DOUBLE, bufB,
ProcPartElem, MPI_DOUBLE, 0, MPI_COMM_WORLD);
temp = 0.0;
for ( i = 0; i
< ProcPartsize; i++)
{
for ( j
= 0; j < ProcPartsize; j++)
{
for
( k = 0; k < dim; k++ )
temp +=
bufA[i*dim+k]*bufB[j*dim+k];
bufC[i*dim+j+ProcPartsize*ProcRank]
= temp;
temp = 0.0;
}
}
int NextProc;
int PrevProc;
for ( p = 1; p
< ProcNum; p++ )
{
NextProc = ProcRank + 1;
if (
ProcRank == ProcNum-1 )
NextProc = 0;
PrevProc = ProcRank - 1;
if (
ProcRank == 0 )
PrevProc =ProcNum-1;
MPI_Sendrecv_replace(bufB,ProcPartElem,
MPI_DOUBLE, NextProc, 0, PrevProc, 0, MPI_COMM_WORLD, &Status);
temp = 0.0;
for ( i
= 0; i < ProcPartsize; i++ )
for
( j = 0; j < ProcPartsize; j++ )
{
for
( k = 0; k < dim; k++ )
temp +=
bufA[i*dim+k]*bufB[j*dim+k];
if
(ProcRank - p >= 0 )
ind=ProcRank - p;
else
ind = ProcNum -
p+ProcRank;
bufC[i*dim+j+ind*ProcPartsize]
= temp;
temp = 0.0;
}
}
MPI_Gather(bufC,
ProcPartElem, MPI_DOUBLE, C, ProcPartElem, MPI_DOUBLE, 0, MPI_COMM_WORLD);
delete []bufA;
delete []bufB;
delete []bufC;
}
Ленточный алгоритм (модификация)
Отличие второго алгоритма состоит в том, что в подзадачах
располагаются не столбцы, а строки матрицы B. Как результат, перемножение
данных каждой подзадачи сводится не к скалярному умножению имеющихся векторов,
а к их поэлементному умножению. В результате подобного умножения в каждой
подзадаче получается строка частичных результатов для матрицы C.
При данном алгоритме разделения данных между потоками для
выполнения операции умножения матриц нужно обеспечить последовательную передачу
данных строк или блоков матрицы B. На каждой итерации данные перемножаются
поэлементно и производится суммирование с полученными ранее значениями.
При рассмотренном способе разделения данных для выполнения
операции матричного умножения нужно обеспечить последовательное получение в
подзадачах всех строк матрицы B, поэлементное умножение данных и суммирование
вновь получаемых значений с ранее вычисленными результатами.
void MatrixMultiplicationMPI2(double *A, double
*B, double *C,int
size)
{
int i = 0, j = 0;
double *bufA,
*bufB, *bufC;
int dim = size;
MPI_Status Status;
int ProcPartsize
= dim/ProcNum;
int ProcPartElem
= ProcPartsize*dim;
bufA = new double[ProcPartElem];
bufB = new double[ProcPartElem];
bufC = new double[ProcPartElem];
for (i = 0; i
< ProcPartElem; i++)
bufC[i] = 0;
MPI_Scatter(A, ProcPartElem, MPI_DOUBLE, bufA,
ProcPartElem, MPI_DOUBLE,
0, MPI_COMM_WORLD);
MPI_Scatter(B, ProcPartElem, MPI_DOUBLE, bufB,
ProcPartElem, MPI_DOUBLE,
0, MPI_COMM_WORLD);
int k = 0 ;
int NextProc =
ProcRank + 1;
if ( ProcRank ==
ProcNum - 1 ) NextProc = 0;
int PrevProc =
ProcRank - 1;
if ( ProcRank ==
0 ) PrevProc = ProcNum - 1;
int PrevDataNum =
ProcRank;
for (int p = 0; p < ProcNum; p++)
{
for (i = 0; i
< ProcPartsize; i++)
for (j = 0;
j < size; j++)
{
double
tmp = 0;
for
(k = 0; k < ProcPartsize; k++)
tmp += bufA[PrevDataNum *
ProcPartsize + i * size + k] * bufB[k * size + j];
bufC[i * size + j] += tmp;
}
PrevDataNum -= 1;
if
(PrevDataNum < 0)
PrevDataNum = ProcNum - 1;
MPI_Sendrecv_replace(bufB, ProcPartElem,
MPI_DOUBLE, NextProc, 0, PrevProc, 0, MPI_COMM_WORLD, &Status);
}
MPI_Gather(bufC, ProcPartElem, MPI_DOUBLE, C,
ProcPartElem, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
Анализ эффективности
Примем за Ts - время работы последовательного
алгоритма, за Tp - время работы параллельного алгоритма, за S -
ускорение, за E - эффективность, за p – количество процессоров.
При применении последовательного алгоритма перемножения
матриц число шагов имеет порядок O(n3).
Для параллельного алгоритма на каждой итерации каждый
процессор выполняет умножение имеющихся на процессоре полос исходных матриц А и
В (размерность каждой полосы составляет n/p, следовательно, общее количество
выполняемых при этом умножении операций равно n3/p2). Поскольку число итераций
алгоритма совпадает с количеством процессоров, сложность параллельного
алгоритма без учета затрат на передачу данных может быть определена при помощи
выражения:
Tp = (n3/p2) * p = n3/p.
Показатели ускорения и эффективности данного параллельного
алгоритма можно записать как:
S = Ts/Tp = n3 / ( n3
/ p) = p
E = n3 / [ ( n3 / p) * p] = 1
Полученные показатели являются идеальными. На практике же мы
обычно получаем относительные показатели, которые с идеальными не совпадают
из-за затрат на выполнение операций передачи данных между процессорами,
служебных затрат и пр.
Оценим реальные показатели с учётом затрат на выполнение
операций передачи данных между процессорами. С учетом числа и длительности
выполняемых операций время выполнения вычислений параллельного алгоритма может
быть оценено следующим образом:
Tp(calc) = [ n2 / p ] · [ 2 n - 1 ] ·
t, где t - время выполнения одной элементарной скалярной операции.
Для оценки коммуникационной сложности параллельных
вычислений будем предполагать, что все операции передачи данных между
процессорами в ходе одной итерации алгоритма могут быть выполнены параллельно.
Объем передаваемых данных между процессорами определяется размером полос и
составляет n/p строк или столбцов длины n. Общее количество параллельных
операций передачи сообщений на единицу меньше числа итераций алгоритма (на
последней итерации передача данных не является обязательной). Тем самым, оценка
трудоемкости выполняемых операций передачи данных может быть определена как:
Tp ( comm ) = ( p - 1 ) · ( a + w · n · ( n / p
) / b ), где a - латентность, b - пропускная способность сети передачи данных,
а w - размер элемента матрицы в байтах. Тогда
Tp = ( n2 / p ) · ( 2 n - 1 ) · t + ( p
- 1 ) · ( a + w · n · ( n / p ) / b ).
Результаты
Эксперимент проводился на следующей системе: Intel Core 2 Duo
2.4 GHz, Microsoft Windows 7 Professional, Компилятор MS Visual Studio 2008.
Размер матрицы |
Последовательный алгоритм |
Параллельный алгоритм 1 (1 процесс) |
Параллельный алгоритм 2 (1 процесс) |
100 |
0.002048 |
0.001446 |
0.001702 |
200 |
0.017264 |
0.010924 |
0.012796 |
300 |
0.063847 |
0.038862 |
0.043494 |
400 |
0.153293 |
0.096983 |
0.282999 |
500 |
0.281513 |
0.23104 |
0.748915 |
600 |
0.498468 |
0.386175 |
1.426799 |
700 |
0.817433 |
0.631975 |
2.477571 |
800 |
1.350454 |
0.981429 |
3.590774 |
900 |
1.852341 |
1.384118 |
5.459218 |
1000 |
2.900795 |
1.979179 |
7.826725 |
Размер матрицы |
Параллельный алгоритм 1 (2 процесса) |
Параллельный алгоритм 2 (2 процесса) |
Ускорение алгоритма 1 |
Ускорение алгоритма 2 |
100 |
0.000933 |
0.001254 |
1.549839228 |
1.357256778 |
200 |
0.008287 |
0.008038 |
1.318209243 |
1.591938293 |
300 |
0.023626 |
0.027977 |
1.644882756 |
1.554634164 |
400 |
0.058177 |
0.100434 |
1.667033364 |
2.817760918 |
500 |
0.118107 |
0.28227 |
1.956192266 |
2.653186665 |
600 |
0.332878 |
0.675594 |
1.16010971 |
2.111917809 |
700 |
0.534014 |
1.219412 |
1.183442756 |
2.031775151 |
800 |
0.914666 |
1.928587 |
1.072991671 |
1.861867782 |
900 |
1.305989 |
3.434378 |
1.059823628 |
1.589579831 |
1000 |
1.651574 |
4.858325 |
1.198359262 |
1.610992472 |
Об авторах
Лазарева Наташа
Тилькова Лена
83-03
