Матрица — математический объект, записываемый в виде прямоугольной таблицы элементов кольца или поля, которая представляет собой совокупность строк и столбцов, на пересечении которых находятся её элементы. Количество строк и столбцов матрицы задают размер матрицы(wiki). Операция умножения матриц является одной из основных задач матричных вычислений. Умножение матрицы A размера heightA*widthA и матрицы B размера heightB*widthB (widthA равно heightB) приводит к получению матрицы С размера heightA*widthB, каждый элемент которой определяется скалярным произведением соответствующей строки матрицы А и столбца матрицы В.

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

1. Последовательный алгоритм

Этот алгоритм предполагает выполнение heightA*widthA*widthB операций умножения и столько же операций сложения элементов исходных матриц. При умножении квадратных матриц размера n*n количество выполненных операций имеет порядок O(n³).

void StandartMatrixMultiplication(int *A, int AMatrixHeght,int AMatrixWidth, int *B, int BMatrixHeght,int BMatrixWidth, int *C)
{
  int tmp;
  for (int i = 0; i < AMatrixHeght; i++)
  {
    for (int j = 0; j < BMatrixWidth; j++)
    {
      C[i * BMatrixWidth + j] = 0; 
      tmp = 0;
      for (int k = 0; k < AMatrixWidth; k ++)
      {
	tmp += A[i * AMatrixWidth + k] * B[k * BMatrixWidth + j];
	C[i * BMatrixWidth + j] = tmp;
      }
    }
  }
}

Этот алгоритм является итеративным и ориентирован на последовательное вычисление строк матрицы С. Действительно, при выполнении одной итерации внешнего цикла (цикла по переменной i) вычисляется одна строка результирующей матрицы.

Поскольку каждый элемент результирующей матрицы есть скалярное произведение строки и столбца исходных матриц, то для вычисления всех элементов матрицы С размером n*n необходимо выполнить n² * ( 2 n - 1 ) скалярных операций и затратить время T1 = n² * ( 2 n - 1 ) * t, где t есть время выполнения одной элементарной скалярной операции.

2. Параллельная схема

Вычисление всех элементов матрицы С = A x B может быть выполнено независимо друг от друга. Исходные матрицы можно разбить на блоки, каждый из которых в свою очередь будет передаваться соответствующему процессу.
Для вычисления строки матрицы С необходимо, чтобы в каждой подзадаче содержалась строка (строки, если передаваемый блок этой матрицы содержит более одной строки или количество процессов как минимум в 2 раза больше чем количество строк в матрице) матрицы А и был обеспечен доступ ко всем столбцам матрицы B.

Возможны различные способы организации параллельных алгоритмов:

Первый алгоритм

Данный алгоритм представляет собой итерационную процедуру. Исходные матрицы разбиваются на блоки. Количество этих блоков зависит от количества потоков в программе. На каждой итерации алгоритма все подзадачи содержат по одной или несколько строк матрицы А и по одному или нескольким столбцам матрицы В. При выполнении итерации проводится скалярное умножение содержащихся в подзадачах блоков друг на друга, что приводит к получению соответствующего блока элементов в вычисляемой матрице С. По завершении вычислений в конце каждой итерации столбцы матрицы В должны быть переданы между подзадачами с тем, чтобы в каждой подзадаче оказались новые блоки этой матрицы для продолжения вычислений матрицы C. На последнем шаге передаются оставшиеся строки матриц A и B. Таким образом, происходит нахождение матрицы C. Количество подзадач зависит от количества строк в блоках передаваемых потокам. Например, если в каждом блоке содержится только по одной строке исходных матриц A и B, то выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью равной n² (по числу элементов матрицы С).

Когда размер матриц n оказывается больше, чем число процессоров p, базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько соседних строк и столбцов перемножаемых матриц. В этом случае исходная матрица A разбивается на ряд горизонтальных полос, а матрица B представляется в виде набора вертикальных полос. Размер полос при этом следует выбрать равным k = n/p (в предположении, что n кратно p), что позволит по-прежнему обеспечить равномерность распределения вычислительной нагрузки по процессорам, составляющим многопроцессорную вычислительную систему.

int  k;
int n_process;
int index;
int tmp;
MPI_Status Status;
int ProcPart_A = Ar_size / ProcNum;
int part_A = ProcPart_A * Ac_size;
int ProcPart_B = Bc_size / ProcNum;
int part_B = ProcPart_B * Ac_size;
int* buf_A = new int[part_A];
int* buf_B = new int[part_B];
int* buf_C = new int[ProcPart_A * Bc_size];
MPI_Scatter(A, part_A, MPI_INT, buf_A, part_A, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Scatter(B, part_B, MPI_INT, buf_B, part_B, MPI_INT, 0, MPI_COMM_WORLD);
tmp = 0;
for(i = 0; i < ProcPart_A; i++)
   for(j = 0; j < ProcPart_B; j++)
   {
      tmp = 0;
      for(k = 0; k < Ac_size; k++)
         tmp += buf_A[i * Ac_size + k] * buf_B[j * Ac_size + k];
      buf_C[i * Bc_size + j + ProcPart_A * ProcRank] = tmp;
   }
int NextProc; 
int PrevProc;
NextProc = ProcRank + 1;
if (ProcRank == ProcNum - 1) NextProc = 0;
PrevProc = ProcRank - 1;
if (ProcRank == 0) PrevProc = ProcNum - 1;
for(n_process = 1; n_process < ProcNum; n_process++) 
{
   MPI_Sendrecv_replace(buf_B, part_B, MPI_INT, NextProc, 0, PrevProc, 0, MPI_COMM_WORLD, &Status);   // передача-принятие единого типа сообщений (рассылка B)
   tmp = 0;
   for (i = 0; i < ProcPart_A; i++)
   {
      for(j = 0; j < ProcPart_B; j++) 
	 {
	    tmp = 0;
            for(k = 0; k < Ac_size; k++)
	       tmp += buf_A[i * Ac_size + k] * buf_B[j * Ac_size + k];     
            if (ProcRank - n_process >= 0)
               index = ProcRank - n_process;
            else 
	       index =(ProcNum - n_process + ProcRank);
            buf_C[i * Bc_size + j + index * ProcPart_A] = tmp;
	 }
   }
}
MPI_Gather(buf_C, ProcPart_A * Bc_size, MPI_INT, C, ProcPart_A * Bc_size, MPI_INT, 0, MPI_COMM_WORLD);

Второй алгоритм

Отличие второго алгоритма состоит в том, что в подзадачах располагаются не столбцы, а строки матрицы B. Как результат, перемножение данных каждой подзадачи сводится не к скалярному умножению имеющихся векторов, а к их поэлементному умножению. В результате подобного умножения в каждой подзадаче получается строка частичных результатов для матрицы C.
При данном алгоритме разделения данных между потоками для выполнения операции умножения матриц нужно обеспечить последовательную передачу данных строк или блоков матрицы B. На каждой итерации данные перемножаются поэлементно и производится суммирование с полученными ранее значениями.

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

void Parallel2MatrixMultiplicftion(int *A, int AMatrixHeght,int AMatrixWidth, int *B, int BMatrixHeght,int BMatrixWidth, int *C, int ProcRank, int ProcNum)
{
   int i = 0, j = 0;
   int *bufA, *bufB, *bufC;
   int rowA = AMatrixHeght / ProcNum;
   int rowB = BMatrixHeght / ProcNum;
   int rowC = rowA;
   int PartA = rowA * AMatrixWidth;
   int PartB = rowB * BMatrixWidth;
   int partC = rowC * BMatrixWidth;
   bufA = new int[PartA];
   bufB = new int[PartB];
   bufC = new int[partC];
   for (i = 0; i < partC; i++) 
      bufC[i] = 0;
   MPI_Scatter(A, PartA, MPI_INT, bufA, PartA, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Scatter(B, PartB, MPI_INT, bufB, PartB, MPI_INT, 0, MPI_COMM_WORLD);
   int k = 0 , tmp = 0;
   int NextProc = ProcRank + 1;
   if ( ProcRank == ProcNum - 1 ) NextProc = 0;
   int PrevProc = ProcRank - 1;
   if ( ProcRank == 0 ) PrevProc = ProcNum - 1;
   MPI_Status Status;
   int PrevDataNum = ProcRank;
   for (int p = 0; p < ProcNum; p++)
   {
      for (i = 0; i < rowA; i++)
      {
         for (j = 0; j < BMatrixWidth; j++)
         {
            tmp = 0;
	    for (k = 0; k < rowB; k++)
	       tmp += bufA[PrevDataNum * rowB + i * AMatrixWidth + k] * bufB[k * BMatrixWidth + j];
	       bufC[i * BMatrixWidth + j] += tmp;
	 }
      }
      PrevDataNum -= 1;
      if (PrevDataNum < 0) PrevDataNum = ProcNum - 1;
      MPI_Sendrecv_replace(bufB, PartB, MPI_INT, NextProc, 0, PrevProc, 0, MPI_COMM_WORLD, &Status);
   }
   MPI_Gather(bufC, partC, MPI_INT, C, partC, MPI_INT, 0, MPI_COMM_WORLD);
}

Примем за T_s - время работы последовательного алгоритма, за T_p - время работы параллельного алгоритма, за S - ускорение, за E - эффективность, за p – количество процессоров.

При применении последовательного алгоритма перемножения матриц число шагов имеет порядок O(n³).
Для параллельного алгоритма на каждой итерации каждый процессор выполняет умножение имеющихся на процессоре полос исходных матриц А и В (размерность каждой полосы составляет n/p, следовательно, общее количество выполняемых при этом умножении операций равно n³/p²). Поскольку число итераций алгоритма совпадает с количеством процессоров, сложность параллельного алгоритма без учета затрат на передачу данных может быть определена при помощи выражения:

Оценим реальные показатели с учётом затрат на выполнение операций передачи данных между процессорами. С учетом числа и длительности выполняемых операций время выполнения вычислений параллельного алгоритма может быть оценено следующим образом:

Для оценки коммуникационной сложности параллельных вычислений будем предполагать, что все операции передачи данных между процессорами в ходе одной итерации алгоритма могут быть выполнены параллельно. Объем передаваемых данных между процессорами определяется размером полос и составляет n/p строк или столбцов длины n. Общее количество параллельных операций передачи сообщений на единицу меньше числа итераций алгоритма (на последней итерации передача данных не является обязательной). Тем самым, оценка трудоемкости выполняемых операций передачи данных может быть определена как:

Первый алгоритм.

Второй алгоритм.

Эксперименты проводились на комплексе из двухпроцессорных вычислительных узлов на базе двухядерных процессоров Intel Core 2 Duo E6550, 2.33 ГГц, 4 Гб RAM, под управлением операционной системы Microsoft Windows 7 Enterprise. Разработка программ проводилась в среде Microsoft Visual Studio 2008, для компиляции использовался стандартный компилятор, предоставляемый средой, с включенной полной оптимизацией.

Умножение матрицы на матрицу при первом ленточном алгоритме разделения данных

Порядок матрицы, элементы	Последовательный алгоритм, секунды	Параллельный алгоритм
		2 процессора		4 процессора		8 процессоров
		Tp, секунды	Sp, секунды	Tp, секунды	Sp, секунды	Tp, секунды	Sp, секунды
64	0,000344	0,000396	0,868686	0,213492	0,001611	0.220349	0,001561
128	0,005160	0,002025	2,548148	0,237700	0,021708	0,230087	0,022426
256	0,041984	0,016802	2,498750	0,271006	0,154919	0,291152	0,144199
512	0,897446	0,128341	6,992667	0,478074	1,877211	1,509614	0,594487
1024	12,640094	1,173534	10,770965	1,453582	8,695824	2,702761	4,676733
2048	221,676410	8,780570	25,246243	8,061305	17,537560	8,623663	25,70559

Умножение матрицы на матрицу при втором ленточном алгоритме разделения данных

Порядок матрицы, элементы	Последовательный алгоритм, секунды	Параллельный алгоритм
		2 процессора		4 процессора		8 процессоров
		Tp, секунды	Sp, секунды	Tp, секунды	Sp, секунды	Tp, секунды	Sp, секунды
64	0,000344	0,000868	0,396313	0,220581	0,001559	0,889506	0,000386
128	0,005160	0,002673	1,930415	0,282324	0,018276	1,236358	0,004173
256	0,041984	0,033341	1,259230	0,932165	0,045039	1,334018	0,031471
512	0,897446	0,269871	3,325462	1,211642	0,740685	2,102328	0,426882
1024	12,640094	3,549540	3,561051	3,698646	3,417492	3,950561	3,199569
2048	221,676410	67,984003	3,2607142	20,244242	10,950096	14,859041	14,918621

Лабораторную работу выполнели студенты группы 8409 факультета ВМК: Панина Катерина и Соломин Сергей