Параллельные методы матричного умножения

7.3.3. Масштабирование и распределение подзадач по процессорам

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

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

7.3.4. Анализ эффективности

Выполним анализ эффективности первого параллельного алгоритма умножения матриц.

Общая трудоемкость последовательного алгоритма, как уже отмечалось ранее, является пропорциональной n3. Для параллельного алгоритма на каждой итерации каждый процессор выполняет умножение имеющихся на процессоре полос матрицы А и матрицы В (размер полос равен n/p, и, как результат, общее количество выполняемых при этом умножении операций равно n3/p2 ). Поскольку число итераций алгоритма совпадает с количеством процессоров, сложность параллельного алгоритма без учета затрат на передачу данных может быть определена при помощи выражения$$T_p = (n^3/p^2)\cdot p = n^3/p .$$

С учетом этой оценки показатели ускорения и эффективности данного параллельного алгоритма матричного умножения принимают вид:$$S_p = \frac{n^3}{n^3/p} = p \quad \text{и} \quad E_p=\frac{n^3}{p\cdot(n^3/p)}=1.$$

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

С учетом числа и длительности выполняемых операций время выполнения вычислений параллельного алгоритма может быть оценено следующим образом:$$T_p(calc)=(n^2/p)\cdot(2n-1)\cdot\tau$$ (здесь, как и ранее, $$\tau$$ есть время выполнения одной элементарной скалярной операции).

Для оценки коммуникационной сложности параллельных вычислений будем предполагать, что все операции передачи данных между процессорами в ходе одной итерации алгоритма могут быть выполнены параллельно. Объем передаваемых данных между процессорами определяется размером полос и составляет n/p строк или столбцов длины n. Общее количество параллельных операций передачи сообщений на единицу меньше числа итераций алгоритма (на последней итерации передача данных не является обязательной). Тем самым, оценка трудоемкости выполняемых операций передачи данных может быть определена как$$T_p(comm)=(p-1)\cdot(\alpha + w\cdot n\cdot(n/p)/\beta),$$ где $$\alpha$$ – латентность, $$\beta$$ – пропускная способность сети передачи данных, а w есть размер элемента матрицы в байтах.

С учетом полученных соотношений общее время выполнения параллельного алгоритма матричного умножения определяется следующим выражением:$$T_p=(n^2/p)(2n-1)\cdot\tau+(p-1)\cdot(\alpha + w\cdot n\cdot(n/p)/\beta).$$

7.3.5. Результаты вычислительных экспериментов

Эксперименты проводились на вычислительном кластере на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server.

Для оценки длительности $$\tau$$ базовой скалярной операции проводилось решение задачи умножения матриц при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины $$\tau$$ было получено значение 6,4 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности a и пропускной способности b соответственно 130 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.

Результаты вычислительных экспериментов приведены в таблице 7.1. Эксперименты выполнялись с использованием двух, четырех и восьми процессоров.

Сравнение экспериментального времени $$T^*_p$$ выполнения эксперимента и теоретического времени Tp из формулы (7.8) представлено в таблице 7.2 и на рис. 7.5.

7.4. Алгоритм Фокса умножения матриц при блочном разделении данных

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

7.4.1. Определение подзадач

Блочная схема разбиения матриц подробно изложена в первом разделе лекции 6. При таком способе разделения данных исходные матрицы А, В и результирующая матрица С представляются в виде наборов блоков. Для более простого изложения следующего материала будем предполагать далее, что все матрицы являются квадратными размера nxn, количество блоков по горизонтали и вертикали одинаково и равно q (т.е. размер всех блоков равен kxk, k=n/q ). При таком представлении данных операция матричного умножения матриц А и B в блочном виде может быть представлена так:$$\begin{pmatrix} A_{00}A_{01}\ldots A_{0q-1} \\ \ldots \\ A_{q-10}A_{q-11}\ldots A_{q-1q-1} \end{pmatrix} \times \begin{pmatrix} B_{00}B_{01}\ldots B_{0q-1} \\ \ldots \\ B_{q-10}B_{q-11}\ldots B_{q-1q-1} \end{pmatrix} = \begin{pmatrix} C_{00}C_{01}\ldots C_{0q-1} \\ \ldots \\ C_{q-10}C_{q-11}\ldots C_{q-1q-1} \end{pmatrix} ,$$ где каждый блок Cij матрицы C определяется в соответствии с выражением$$C_{ij}=\sum_{s=0}^{q-1} A_{is} B_{sj}.$$

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

Для выполнения всех необходимых вычислений базовым подзадачам должны быть доступны соответствующие наборы строк матрицы A и столбцов матрицы B. Размещение всех требуемых данных в каждой подзадаче неизбежно приведет к дублированию и к значительному росту объема используемой памяти. Как результат, вычисления должны быть организованы таким образом, чтобы в каждый текущий момент времени подзадачи содержали лишь часть необходимых для проведения расчетов данных, а доступ к остальной части данных обеспечивался бы при помощи передачи данных между процессорами. Один из возможных подходов – алгоритм Фокса ( Fox ) – рассмотрен далее в данном подразделе. Второй способ – алгоритм Кэннона ( Cannon ) – приводится в подразделе 7.5.

7.4.2. Выделение информационных зависимостей

Итак, за основу параллельных вычислений для матричного умножения при блочном разделении данных принят подход, при котором базовые подзадачи отвечают за вычисления отдельных блоков матрицы C и при этом в подзадачах на каждой итерации расчетов располагается только по одному блоку исходных матриц A и B. Для нумерации подзадач будем использовать индексы размещаемых в подзадачах блоков матрицы C, т.е. подзадача (i,j) отвечает за вычисление блока Cij – тем самым, набор подзадач образует квадратную решетку, соответствующую структуре блочного представления матрицы C.

Возможный способ организации вычислений при таких условиях состоит в применении широко известного алгоритма Фокса ( Fox ) — см., например, [, ].

В соответствии с алгоритмом Фокса в ходе вычислений на каждой базовой подзадаче (i,j) располагается четыре матричных блока:

Выполнение параллельного метода включает:

Для пояснения этих правил параллельного метода на рис. 7.6 приведено состояние блоков в каждой подзадаче в ходе выполнения итераций этапа вычислений (для решетки подзадач 2x2 ).

7.4.3. Масштабирование и распределение подзадач по процессорам

В рассмотренной схеме параллельных вычислений количество блоков может варьироваться в зависимости от выбора размера блоков – эти размеры могут быть подобраны таким образом, чтобы общее количество базовых подзадач совпадало с числом процессоров p. Так, например, в наиболее простом случае, когда число процессоров представимо в виде $$p=\delta 2$$ (т.е. является полным квадратом), можно выбрать количество блоков в матрицах по вертикали и горизонтали равным $$\delta$$ (т.е. $$q=\delta$$ ). Такой способ определения количества блоков приводит к тому, что объем вычислений в каждой подзадаче является одинаковым и тем самым достигается полная балансировка вычислительной нагрузки между процессорами. В более общем случае при произвольных количестве процессоров и размерах матриц балансировка вычислений может отличаться от абсолютно одинаковой, но, тем не менее, при надлежащем выборе параметров может быть распределена между процессорами равномерно в рамках требуемой точности.

Для эффективного выполнения алгоритма Фокса, в котором базовые подзадачи представлены в виде квадратной решетки и в ходе вычислений выполняются операции передачи блоков по строкам и столбцам решетки подзадач, наиболее адекватным решением является организация множества имеющихся процессоров также в виде квадратной решетки. В этом случае можно осуществить непосредственное отображение набора подзадач на множество процессоров – базовую подзадачу (i,j) следует располагать на процессоре Pi,j. Необходимая структура сети передачи данных может быть обеспечена на физическом уровне, если топология вычислительной системы имеет вид решетки или полного графа.

7.4.4. Анализ эффективности

Определим вычислительную сложность данного алгоритма Фокса. Построение оценок будет происходить при условии выполнения всех ранее выдвинутых предположений: все матрицы являются квадратными размера nxn, количество блоков по горизонтали и вертикали являются одинаковым и равным q (т.е. размер всех блоков равен kxk, k=n/q ), процессоры образуют квадратную решетку и их количество равно p=q2.

Как уже отмечалось, алгоритм Фокса требует для своего выполнения q итераций, в ходе которых каждый процессор перемножает свои текущие блоки матриц А и В и прибавляет результаты умножения к текущему значению блока матрицы C. С учетом выдвинутых предположений общее количество выполняемых при этом операций будет иметь порядок n3/p. Как результат, показатели ускорения и эффективности алгоритма имеют вид:$$S_p = \frac{n^3}{n^3/p} = p \quad \text{и} \quad E_p=\frac{n^3}{p\cdot(n^3/p)}=1.$$

Общий анализ сложности снова дает идеальные показатели эффективности параллельных вычислений. Уточним полученные соотношения — для этого укажем более точно количество вычислительных операций алгоритма и учтем затраты на выполнение операций передачи данных между процессорами.

Определим количество вычислительных операций. Сложность выполнения скалярного умножения строки блока матрицы A на столбец блока матрицы В можно оценить как 2(n/q)-1. Количество строк и столбцов в блоках равно n/q и, как результат, трудоемкость операции блочного умножения оказывается равной (n2/p)(2n/q-1). Для сложения блоков требуется n2/p операций. С учетом всех перечисленных выражений время выполнения вычислительных операций алгоритма Фокса может быть оценено следующим образом:$$T_p(calc)=q[(n^2 /p)\cdot(2n/q-1)+(n^2 /p)]\cdot\tau.$$ (напомним, что $$\tau$$ есть время выполнения одной элементарной скалярной операции).

Оценим затраты на выполнение операций передачи данных между процессорами. На каждой итерации алгоритма перед умножением блоков один из процессоров строки процессорной решетки рассылает свой блок матрицы A остальным процессорам своей строки. Как уже отмечалось ранее, при топологии сети в виде гиперкуба или полного графа выполнение этой операции может быть обеспечено за log2q шагов, а объем передаваемых блоков равен n2/p. Как результат, время выполнения операции передачи блоков матрицы A при использовании модели Хокни может оцениваться как$$T_p^1(comm)=\log_2 q(\alpha +w(n^2/p)/ \beta),$$ где $$\alpha$$ – латентность, $$\beta$$ – пропускная способность сети передачи данных, а w есть размер элемента матрицы в байтах. В случае же когда топология строк процессорной решетки представляет собой кольцо, выражение для оценки времени передачи блоков матрицы A принимает вид:$$\widetilde{T}_p^1(comm)=(q/2)(\alpha+w(n^2/p)/ \beta).$$

Далее после умножения матричных блоков процессоры передают свои блоки матрицы В предыдущим процессорам по столбцам процессорной решетки (первые процессоры столбцов передают свои данные последним процессорам в столбцах решетки). Эти операции могут быть выполнены процессорами параллельно, и, тем самым, длительность такой коммуникационной операции составляет:$$T-p^2(comm)=\alpha+w \cdot (n^2/p)/ \beta.$$

Просуммировав все полученные выражения, можно получить, что общее время выполнения алгоритма Фокса может быть определено при помощи следующих соотношений:$$\begin{aligned} T_p = q[(n^2/p) \cdot (2n/q-1)+(n^2/p)] \cdot \tau+q \log_2 q(\alpha + w(n^2/p)/ \beta) + \\ +(q-1) \cdot (\alpha + w(n^2/p)/ \beta) = q[(n^2/p) \cdot (2n/q-1)+(n^2/p)] \cdot \tau+ \\ +(q\log_2 q +(q-1))( \alpha + w(n^2/p)/ \beta) \end{aligned}$$ (напомним, что параметр q определяет размер процессорной решетки и $$q=\sqrt{p}$$ ).

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

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

1. Главная функция программы. Определяет основную логику работы алгоритма, последовательно вызывает необходимые подпрограммы.

// Программа 7.1
// Алгоритм Фокса умножения матриц – блочное представление данных
//   Условия выполнения программы: все матрицы квадратные, 
//   размер блоков и их количество по горизонтали и вертикали
//   одинаково, процессы образуют квадратную решетку
int ProcNum = 0;   	// Количество доступных процессов 
int ProcRank = 0;  	// Ранг текущего процесса
int GridSize;      	// Размер виртуальной решетки процессов
int GridCoords[2];  	// Координаты текущего процесса в процессной
// решетке
MPI_Comm GridComm;  	// Коммуникатор в виде квадратной решетки
MPI_Comm ColComm;   	// коммуникатор – столбец решетки
MPI_Comm RowComm;   	// коммуникатор – строка решетки

void main ( int argc, char * argv[] ) {
  double* pAMatrix; 	// Первый аргумент матричного умножения
  double* pBMatrix; 	// Второй аргумент матричного умножения
  double* pCMatrix; 	// Результирующая матрица
  int Size;        	// Размер матриц
  int BlockSize;   	// Размер матричных блоков, расположенных 
                  	// на процессах
  double *pAblock;  	// Блок матрицы А на процессе
  double *pBblock;  	// Блок матрицы В на процессе
  double *pCblock;  	// Блок результирующей матрицы С на процессе
  double *pMatrixAblock;
  double Start, Finish, Duration;

  setvbuf(stdout, 0, _IONBF, 0);

  MPI_Init(argc, argv);
  MPI_Comm_size(MPI_COMM_WORLD, ProcNum);
  MPI_Comm_rank(MPI_COMM_WORLD, ProcRank);

  GridSize = sqrt((double)ProcNum);
  if (ProcNum != GridSize*GridSize) {
    if (ProcRank == 0) {
      printf ("Number of processes must be a perfect square \n");
    }
  }
  else {
    if (ProcRank == 0)
      printf("Parallel matrix multiplication program\n");

    // Создание виртуальной решетки процессов и коммуникаторов 
    // строк и столбцов
    CreateGridCommunicators();
  
    // Выделение памяти и инициализация элементов матриц
    ProcessInitialization ( pAMatrix, pBMatrix, pCMatrix, pAblock,
      pBblock, pCblock, pMatrixAblock, Size, BlockSize );
    // Блочное распределение матриц между процессами
    DataDistribution(pAMatrix, pBMatrix, pMatrixAblock, pBblock, Size, 
      BlockSize);

    // Выполнение параллельного метода Фокса
    ParallelResultCalculation(pAblock, pMatrixAblock, pBblock, 
      pCblock, BlockSize);

    // Сбор результирующей матрицы на ведущем процессе
    ResultCollection(pCMatrix, pCblock, Size, BlockSize);

    // Завершение процесса вычислений
    ProcessTermination (pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock, 
      pCblock, pMatrixAblock);
  }

  MPI_Finalize();
}

2. Функция CreateGridCommunicators. Данная функция создает коммуникатор в виде двумерной квадратной решетки, определяет координаты каждого процесса в этой решетке, а также создает коммуникаторы отдельно для каждой строки и каждого столбца.

Создание решетки производится при помощи функции MPI_Cart_create (вектор Periodic определяет возможность передачи сообщений между граничными процессами строк и столбцов создаваемой решетки). После создания решетки каждый процесс параллельной программы будет иметь координаты своего положения в решетке; получение этих координат обеспечивается при помощи функции MPI_Cart_coords.

Формирование топологий завершается созданием множества коммуникаторов для каждой строки и каждого столбца решетки в отдельности (функция MPI_Cart_sub ).

// Создание коммуникатора в виде двумерной квадратной решетки 
// и коммуникаторов для каждой строки и каждого столбца решетки
void CreateGridCommunicators() {
  int DimSize[2];  	// Количество процессов в каждом измерении 
                  	// решетки
  int Periodic[2];	// =1 для каждого измерения, являющегося 
                 	// периодическим 
  int Subdims[2];  	// =1 для каждого измерения, оставляемого 
                  	// в подрешетке
  DimSize[0] = GridSize; 
  DimSize[1] = GridSize;
  Periodic[0] = 0;
  Periodic[1] = 0;

  // Создание коммуникатора в виде квадратной решетки 
  MPI_Cart_create(MPI_COMM_WORLD, 2, DimSize, Periodic, 1, GridComm);

  // Определение координат процесса в решетке 
  MPI_Cart_coords(GridComm, ProcRank, 2, GridCoords);
  
  // Создание коммуникаторов для строк процессной решетки
  Subdims[0] = 0;  // Фиксация измерения
  Subdims[1] = 1;  // Наличие данного измерения в подрешетке
  MPI_Cart_sub(GridComm, Subdims, RowComm);
  
  // Создание коммуникаторов для столбцов процессной решетки
  Subdims[0] = 1;
  Subdims[1] = 0;
  MPI_Cart_sub(GridComm, Subdims, ColComm);
}

3. Функция ProcessInitialization. Данная функция определяет параметры решаемой задачи (размеры матриц и их блоков), выделяет память для хранения данных и осуществляет ввод исходных матриц (или формирует их при помощи какого-либо датчика случайных чисел). Всего в каждом процессе должна быть выделена память для хранения четырех блоков – для указателей на выделенную память используются переменные pAblock, pBblock, pCblock, pMatrixAblock. Первые три указателя определяют блоки матриц A, B и C соответственно. Следует отметить, что содержимое блоков pAblock и pBblock постоянно меняется в соответствии с пересылкой данных между процессами, в то время как блок pMatrixAblock матрицы A остается неизменным и применяется при рассылках блоков по строкам решетки процессов (см. функцию AblockCommunication ).

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

// Функция для выделения памяти и инициализации исходных данных
void ProcessInitialization (double* pAMatrix, double* pBMatrix, 
  double* pCMatrix, double* pAblock, double* pBblock, 
  double* pCblock, double* pTemporaryAblock, int Size, 
  int BlockSize ) {
  if (ProcRank == 0) {
    do {
      printf("\nВведите размер матриц: ");
      scanf("%d", Size);
  
      if (Size%GridSize != 0) {
        printf ("Размер матриц должен быть кратен размеру сетки! \n");
      }
    }
    while (Size%GridSize != 0);
  }
  MPI_Bcast(Size, 1, MPI_INT, 0, MPI_COMM_WORLD);

  BlockSize = Size/GridSize;

  pAblock = new double [BlockSize*BlockSize];
  pBblock = new double [BlockSize*BlockSize];
  pCblock = new double [BlockSize*BlockSize];
  pTemporaryAblock = new double [BlockSize*BlockSize];

  for (int i=0; i<BlockSize*BlockSize; i++) {
    pCblock[i] = 0;
  }
  if (ProcRank == 0) {
    pAMatrix = new double [Size*Size];
    pBMatrix = new double [Size*Size];
    pCMatrix = new double [Size*Size];
    RandomDataInitialization(pAMatrix, pBMatrix, Size);
  } 
}

4. Функции DataDistribution и ResultCollection. После задания исходных матриц на нулевом процессе необходимо осуществить распределение исходных данных. Для этого предназначена функция DataDistribution. Может быть предложено два способа выполнения блочного разделения матриц между процессорами, организованными в двумерную квадратную решетку. В первом из них для организации передачи блоков в рамках одной и той же коммуникационной операции можно сформировать средствами MPI производный тип данных. Во втором способе можно организовать двухэтапную процедуру. На первом этапе матрица разделяется на горизонтальные полосы. Эти полосы распределяются на процессы, составляющие нулевой столбец процессорной решетки. Далее каждая полоса разделяется на блоки между процессами, составляющими строки процессорной решетки.

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

Реализация функций DataDistribution и ResultCollection представляет собой задание для самостоятельной работы.

5. Функция AblockCommunication. Функция выполняет рассылку блоков матрицы A по строкам процессорной решетки. Для этого в каждой строке решетки определяется ведущий процесс Pivot, осуществляющий рассылку. Для рассылки используется блок pMatrixAblock, переданный в процесс в момент начального распределения данных. Выполнение операции рассылки блоков осуществляется при помощи функции MPI_Bcast. Следует отметить, что данная операция является коллективной и ее локализация пределами отдельных строк решетки обеспечивается за счет использования коммуникаторов RowComm, определенных для набора процессов каждой строки решетки в отдельности.

// Рассылка блоков матрицы А по строкам решетки процессов 
void ABlockCommunication (int iter, double *pAblock, 
  double* pMatrixAblock, int BlockSize) {

  // Определение ведущего процесса в строке процессной решетки 
  int Pivot = (GridCoords[0] + iter) % GridSize;
  
  // Копирование передаваемого блока в отдельный буфер памяти
  if (GridCoords[1] == Pivot) {
    for (int i=0; i<BlockSize*BlockSize; i++)
      pAblock[i] = pMatrixAblock[i];
  }
  
  // Рассылка блока
  MPI_Bcast(pAblock, BlockSize*BlockSize, MPI_DOUBLE, Pivot,
    RowComm);
}

6. Функция BlockMultiplication. Функция обеспечивает перемножение блоков матриц A и B. Следует отметить, что для более легкого понимания рассматриваемой программы приводится простой вариант реализации функции – выполнение операции блочного умножения может быть существенным образом оптимизировано для сокращения времени вычислений. Данная оптимизация может быть направлена, например, на повышение эффективности использования кэша процессоров, векторизации выполняемых операций и т.п.

// Умножение матричных блоков
void BlockMultiplication (double *pAblock, double *pBblock, 
  double *pCblock, int BlockSize) {
  // Вычисление произведения матричных блоков
  for (int i=0; i<BlockSize; i++) {
    for (int j=0; j<BlockSize; j++) {
      double temp = 0;
      for (int k=0; k<BlockSize; k++ )
        temp += pAblock [i*BlockSize + k] * pBblock [k*BlockSize + j];
      pCblock [i*BlockSize + j] += temp;
    }
  }
}

7. Функция BblockCommunication. Функция выполняет циклический сдвиг блоков матрицы B по столбцам процессорной решетки. Каждый процесс передает свой блок следующему процессу NextProc в столбце процессов и получает блок, переданный из предыдущего процесса PrevProc в столбце решетки. Выполнение операций передачи данных осуществляется при помощи функции MPI_SendRecv_replace, которая обеспечивает все необходимые пересылки блоков, используя при этом один и тот же буфер памяти pBblock. Кроме того, эта функция гарантирует отсутствие возможных тупиков, когда операции передачи данных начинают одновременно выполняться несколькими процессами при кольцевой топологии сети.

// Циклический сдвиг блоков матрицы В вдоль столбца процессной 
// решетки 
void BblockCommunication (double *pBblock, int BlockSize) {
  MPI_Status Status;
  int NextProc = GridCoords[0] + 1;
  if ( GridCoords[0] == GridSize-1 ) NextProc = 0;
  int PrevProc = GridCoords[0] - 1;
  if ( GridCoords[0] == 0 ) PrevProc = GridSize-1;
  MPI_Sendrecv_replace( pBblock, BlockSize*BlockSize, MPI_DOUBLE,
    NextProc, 0, PrevProc, 0, ColComm, Status);
}

8. Функция ParallelResultCalculation. Для непосредственного выполнения параллельного алгоритма Фокса умножения матриц предназначена функция ParallelResultCalculation, которая реализует логику работы алгоритма.

// Функция для параллельного умножения матриц
void ParallelResultCalculation(double* pAblock, double* pMatrixAblock, 
   double* pBblock, double* pCblock, int BlockSize) {
  for (int iter = 0; iter < GridSize; iter ++) {
    // Рассылка блоков матрицы A по строкам процессной решетки
    ABlockCommunication (iter, pAblock, pMatrixAblock, BlockSize);
    // Умножение блоков
    BlockMultiplication(pAblock, pBblock, pCblock, BlockSize);
    // Циклический сдвиг блоков матрицы B в столбцах процессной 
    // решетки
    BblockCommunication(pBblock, BlockSize);
  }
}

7.4.6. Результаты вычислительных экспериментов

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

Сравнение времени выполнения эксперимента $$T^*_p$$ и теоретического времени Tp, вычисленного в соответствии с выражением (7.13), представлено в таблице 7.4 и на рис. 7.8.

7.5. Алгоритм Кэннона умножения матриц при блочном разделении данных

Рассмотрим еще один параллельный алгоритм матричного умножения, основанный на блочном разбиении матриц, – алгоритм Кэннона ( Cannon ).

7.5.1. Определение подзадач

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

7.5.2. Выделение информационных зависимостей

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

С учетом высказанных замечаний этап инициализации алгоритма Кэннона включает выполнение следующих операций передач данных:

Выполняемые при перераспределении матричных блоков процедуры передачи данных являются примером операции циклического сдвига – см. лекцию 3. Для пояснения используемого способа начального распределения данных на рис. 7.9 показан пример расположения блоков для решетки подзадач 3і3.

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

7.5.3. Масштабирование и распределение подзадач по процессорам

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

Для распределения подзадач между процессорами может быть применен подход, использованный в алгоритме Фокса, – множество имеющихся процессоров представляется в виде квадратной решетки и размещение базовых подзадач (i,j) осуществляется на процессорах Pi,j соответствующих узлов процессорной решетки. Необходимая структура сети передачи данных, как и ранее, может быть обеспечена на физическом уровне при топологии вычислительной системы в виде решетки или полного графа.

7.5.4. Анализ эффективности

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

В соответствии с правилами алгоритма на этапе инициализации производится перераспределение блоков матриц А и B при помощи циклического сдвига матричных блоков по строкам и столбцам процессорной решетки. Трудоемкость выполнения такой операции передачи данных существенным образом зависит от топологии сети. Для сети со структурой полного графа все необходимые пересылки блоков могут быть выполнены одновременно (т.е. длительность операции оказывается равной времени передачи одного матричного блока между соседними процессорами). Для сети с топологией гиперкуба операция циклического сдвига может потребовать выполнения log2q итераций. Для сети с кольцевой структурой связей необходимое количество итераций оказывается равным q-1 – более подробно методы выполнения операции циклического сдвига рассмотрены в лекции 3. Используем для построения оценки коммуникационной сложности этапа инициализации вариант топологии полного графа как более соответствующего кластерным вычислительным системам, время выполнения начального перераспределения блоков может оцениваться как$$T_p^1(comm)=2\cdot(\alpha+w\cdot(n^2/p)/ \beta)$$ (выражение n2/p определяет размер пересылаемых блоков, а коэффициент 2 соответствует двум выполняемым операциям циклического сдвига).

Оценим теперь затраты на передачу данных между процессорами при выполнении основной части алгоритма Кэннона. На каждой итерации алгоритма после умножения матричных блоков процессоры передают свои блоки предыдущим процессорам по строкам (для блоков матрицы A ) и столбцам (для блоков матрицы B ) процессорной решетки. Эти операции также могут быть выполнены процессорами параллельно, и, тем самым, длительность таких коммуникационных действий составляет:$$T_p^2(comm)=2\cdot(\alpha+w\cdot(n^2/p)/ \beta)$$

Поскольку количество итераций алгоритма Кэннона равно q, то с учетом оценки (7.10) общее время выполнения параллельных вычислений может быть определено при помощи следующего соотношения:$$T_p=q[(n^2/p)\cdot(2n/q-1)+(n^2/p)]\cdot\tau+(2q+2)(\alpha+w(n^2/p)/ \beta)$$ (в используемых выражениях параметр $$q=\sqrt{p}$$ определяет размер процессорной решетки).

7.5.5. Результаты вычислительных экспериментов

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

Сравнение времени выполнения эксперимента $$T^*_p$$ и теоретического времени Tp, вычисленного в соответствии с выражением (7.16), представлено в таблице 7.6 и на рис. 7.11.

7.6. Краткий обзор лекции

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

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

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

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

7.7. Обзор литературы

Задача умножения матриц широко рассматривается в литературе. В качестве дополнительного учебного материала могут быть рекомендованы работы [, , ]. Широкое обсуждение вопросов параллельного выполнения матричных вычислений приводится в работе [].

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

7.8. Контрольные вопросы

7.9. Задачи и упражнения