Параллельные методы умножения матрицы на вектор

6.2. Постановка задачи

В результате умножения матрицы А размерности m x n и вектора b, состоящего из n элементов, получается вектор c размера m, каждый i -й элемент которого есть результат скалярного умножения i -й строки матрицы А (обозначим эту строчку ai ) и вектора b:$$c_i = (a_i, b) = \sum_{j=1}^n a_{ij}b_j, 1 \le i \le m.$$

Тем самым получение результирующего вектора c предполагает повторение m однотипных операций по умножению строк матрицы A и вектора b. Каждая такая операция включает умножение элементов строки матрицы и вектора b ( n операций) и последующее суммирование полученных произведений ( n-1 операций). Общее количество необходимых скалярных операций есть величина

T1=mx(2n-1)

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

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

Алгоритм 6.1. Последовательный алгоритм умножения матрицы на вектор

// Алгоритм 6.1
// Поcледовательный алгоритм умножения матрицы на вектор
for (i = 0; i < m; i++){
  c[i] = 0;
  for (j = 0; j < n; j++){
    c[i] += A[i][j]*b[j]
  }
}

Матрично-векторное умножение – это последовательность вычисления скалярных произведений. Поскольку каждое вычисление скалярного произведения векторов длины n требует выполнения n операций умножения и n-1 операций сложения, его трудоемкость порядка O(n). Для выполнения матрично-векторного умножения необходимо осуществить m операций вычисления скалярного произведения, таким образом, алгоритм имеет трудоемкость порядка O(mn).

6.4. Разделение данных

При выполнении параллельных алгоритмов умножения матрицы на вектор, кроме матрицы А, необходимо разделить еще вектор b и вектор результата c. Элементы векторов можно продублировать, то есть скопировать все элементы вектора на все процессоры, составляющие многопроцессорную вычислительную систему, или разделить между процессорами. При блочном разбиении вектора из n элементов каждый процессор обрабатывает непрерывную последовательность из k элементов вектора (мы предполагаем, что размерность вектора n нацело делится на число процессоров, т.е. n= kxp ).

Поясним, почему дублирование векторов b и c между процессорами является допустимым решением (далее для простоты изложения будем полагать, что m=n ). Векторы b и с состоят из n элементов, т.е. содержат столько же данных, сколько и одна строка или один столбец матрицы. Если процессор хранит строку или столбец матрицы и одиночные элементы векторов b и с, то общее число сохраняемых элементов имеет порядок O(n). Если процессор хранит строку (столбец) матрицы и все элементы векторов b и с, то общее число сохраняемых элементов также порядка O(n). Таким образом, при дублировании и при разделении векторов требования к объему памяти из одного класса сложности.

6.5. Умножение матрицы на вектор при разделении данных по строкам

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

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

Для выполнения базовой подзадачи скалярного произведения процессор должен содержать соответствующую строку матрицы А и копию вектора b. После завершения вычислений каждая базовая подзадача определяет один из элементов вектора результата c.

Для объединения результатов расчета и получения полного вектора c на каждом из процессоров вычислительной системы необходимо выполнить операцию обобщенного сбора данных (см. лекцию 4), в которой каждый процессор передает свой вычисленный элемент вектора c всем остальным процессорам. Этот шаг можно выполнить, например, с использованием функции MPI_Allgather из библиотеки MPI.

В общем виде схема информационного взаимодействия подзадач в ходе выполняемых вычислений показана на рис. 6.2.

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

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

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

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

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

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

Рассмотрим трудоемкость алгоритма умножения матрицы на вектор. В случае если матрица А квадратная ( m=n ), последовательный алгоритм умножения матрицы на вектор имеет сложность T1=n2. В случае параллельных вычислений каждый процессор производит умножение только части (полосы) матрицы A на вектор b, размер этих полос равен n/p строк. При вычислении скалярного произведения одной строки матрицы и вектора необходимо произвести n операций умножения и (n-1) операций сложения. Следовательно, вычислительная трудоемкость параллельного алгоритма определяется выражением:$$T_p=n^2/p.$$

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

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

(здесь и далее операция $$\lcei\ \rceil$$ есть округление до целого в большую сторону).

Оценка трудоемкости операции обобщенного сбора данных уже выполнялась в лекции 4 (см. п. 4.3.4). Как уже отмечалась ранее, данная операция может быть выполнена за $$\lceil$$ log2p $$\rceil$$ итерацийБудем предполагать, что топология вычислительной системы допускает возможность такого эффективного способа выполнения операции обобщенного сбора данных (это возможно, в частности, если структура сети передачи данных имеет вид гиперкуба или полного графа).. На первой итерации взаимодействующие пары процессоров обмениваются сообщениями объемом $$w \lceil n/p\rceil$$ ( w есть размер одного элемента вектора c в байтах), на второй итерации этот объем увеличивается вдвое и оказывается равным $$2w\lceil n/p \rceil$$ и т.д. Как результат, длительность выполнения операции сбора данных при использовании модели Хокни может быть определена при помощи следующего выражения$$T_p(comm) = \sum_{i=1}^{\lceil \log_2 p \rceil} (\alpha + 2^{i-1}w\lceil n/p \rceil / \beta) = \alpha \lceil \log_2 p \rceil + w \lceil n/p \rceil (2^{\lceil \log_2 p \rceil} - 1)/\beta,$$ где $$\alpha$$ – латентность сети передачи данных, $$\beta$$ – пропускная способность сети. Таким образом, общее время выполнения параллельного алгоритма составляет$$T_p=(n/p)\cdot (2n-1) \cdot \tau + \alpha \cdot \log_2 p +w(n/p)(p-1)/ \beta$$ (для упрощения выражения в (6.8) предполагалось, что значения n/p и log2p являются целыми).

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

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

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

// Программа 6.1
// Умножение матрицы на вектор – ленточное горизонтальное разбиение
// (исходный и результирующий векторы дублируются между процессами)
void main(int argc, char* argv[]) {
  double* pMatrix;  // Первый аргумент – исходная матрица
  double* pVector;  // Второй аргумент – исходный вектор
  double* pResult;  // Результат умножения матрицы на вектор 
  int Size;        // Размеры исходных матрицы и вектора
  double* pProcRows;
  double* pProcResult;
  int RowNum;
  double Start, Finish, Duration;

  MPI_Init(argc, argv);
  MPI_Comm_size(MPI_COMM_WORLD, ProcNum);
  MPI_Comm_rank(MPI_COMM_WORLD, ProcRank);
  
  // Выделение памяти и инициализация исходных данных
  ProcessInitialization(pMatrix, pVector, pResult, pProcRows, 
                        pProcResult, Size, RowNum);

  // Распределение исходных данных между процессами
  DataDistribution(pMatrix, pProcRows, pVector, Size, RowNum);

  // Параллельное выполнение умножения матрицы на вектор
  ParallelResultCalculation(pProcRows, pVector, pProcResult, 
                            Size, RowNum);

  // Сбор результирующего вектора на всех процессах
  ResultReplication(pProcResult, pResult, Size, RowNum);

  // Завершение процесса вычислений
  ProcessTermination(pMatrix, pVector, pResult, pProcRows, 
                     pProcResult);

  MPI_Finalize();
}

2. Функция ProcessInitialization. Эта функция задает размер и элементы для матрицы A и вектора b. Значения для матрицы A и вектора b определяются в функции RandomDataInitialization.

// Функция для выделения памяти и инициализации исходных данных
void ProcessInitialization (double* pMatrix, double* pVector, 
  double* pResult, double* pProcRows, double* pProcResult, 
  int Size, int RowNum) {
  int RestRows;	// Количество строк матрицы, которые еще 
                // не распределены
  int i;

  if (ProcRank == 0) {
    do {
      printf("\nВведите размер матрицы: ");
      scanf("%d", Size);
      if (Size < ProcNum) {
        printf("Размер матрицы должен превышать количество 
процессов! \n ");
      }
    }
    while (Size < ProcNum);
  }
  MPI_Bcast(Size, 1, MPI_INT, 0, MPI_COMM_WORLD);

  RestRows = Size;
  for (i=0; i<ProcRank; i++) 
    RestRows = RestRows-RestRows/(ProcNum-i);
  RowNum = RestRows/(ProcNum-ProcRank);

  pVector = new double [Size];
  pResult = new double [Size];
  pProcRows = new double [RowNum*Size];
  pProcResult = new double [RowNum];

  if (ProcRank == 0) {
    pMatrix = new double [Size*Size];
    RandomDataInitialization(pMatrix, pVector, Size);
  }
}

3. Функция DataDistribution. Осуществляет рассылку вектора b и распределение строк исходной матрицы A по процессам вычислительной системы. Следует отметить, что когда количество строк матрицы n не является кратным числу процессоров p, объем пересылаемых данных для процессов может оказаться разным и для передачи сообщений необходимо использовать функцию MPI_Scatterv библиотеки MPI.

// Функция для распределения исходных данных между процессами
void DataDistribution(double* pMatrix, double* pProcRows, 
  double* pVector, int Size, int RowNum) {
  int *pSendNum;    	// Количество элементов, посылаемых процессу
  int *pSendInd;    	// Индекс первого элемента данных, 
			// посылаемого процессу
  int RestRows=Size;	// Количество строк матрицы, которые еще 
			// не распределены

  MPI_Bcast(pVector, Size, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  // Выделение памяти для хранения временных объектов
  pSendInd = new int [ProcNum];
  pSendNum = new int [ProcNum];

  // Определение положения строк матрицы, предназначенных 
  // каждому процессу
  RowNum = (Size/ProcNum);
  pSendNum[0] = RowNum*Size;
  pSendInd[0] = 0;
  for (int i=1; i<ProcNum; i++) {
    RestRows -= RowNum;
    RowNum = RestRows/(ProcNum-i);
    pSendNum[i] = RowNum*Size;
    pSendInd[i] = pSendInd[i-1]+pSendNum[i-1];
  }
  // Рассылка строк матрицы
  MPI_Scatterv(pMatrix, pSendNum, pSendInd, MPI_DOUBLE, pProcRows, 
    pSendNum[ProcRank], MPI_DOUBLE, 0, MPI_COMM_WORLD);

  // Освобождение памяти
  delete [] pSendNum;
  delete [] pSendInd;                
}

Следует отметить, что такое разделение действий генерации исходных данных и их рассылки между процессами может быть неоправданным в реальных параллельных вычислениях при большом объеме данных. Широко используемый подход в таких случаях состоит в организации передачи сообщений процессам сразу же после того, как данные процессов будут сгенерированы. Снижение затрат памяти для хранения данных может быть достигнуто также и за счет организации генерации данных в последнем процессе (при таком подходе память для пересылаемых данных и для данных процесса может быть одной и той же).

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

// Функция для вычисления части результирующего вектора
void ParallelResultCalculation(double* pProcRows, double* pVector, 
  double* pProcResult, int Size, int RowNum) {
  int i, j;
  for (i=0; i<RowNum; i++) {
    pProcResult[i] = 0;
    for (j=0; j<Size; j++)
      pProcResult[i] += pProcRows[i*Size+j]*pVector[j];
  }
}

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

// Функция для сбора результирующего вектора на всех процессах
void ResultReplication(double* pProcResult, double* pResult, 
int Size, int RowNum) {
  int *pReceiveNum;  // Количество элементов, посылаемых процессом
  int *pReceiveInd;  // Индекс элемента данных в результирующем 
                     // векторе
  int RestRows=Size; // Количество строк матрицы, которые еще не 
                     // распределены
  int i;

  // Выделение памяти для временных объектов
  pReceiveNum = new int [ProcNum];
  pReceiveInd = new int [ProcNum];

  // Определение положения блоков результирующего вектора 
  pReceiveInd[0] = 0;
  pReceiveNum[0] = Size/ProcNum;
  for (i=1; i<ProcNum; i++) {
    RestRows -= pReceiveNum[i-1];
    pReceiveNum[i] = RestRows/(ProcNum-i);
    pReceiveInd[i] = pReceiveInd[i-1]+pReceiveNum[i-1];
  }
  // Сбор всего результирующего вектора на всех процессах
  MPI_Allgatherv(pProcResult, pReceiveNum[ProcRank], 
    MPI_DOUBLE, pResult, pReceiveNum, pReceiveInd, 
    MPI_DOUBLE, MPI_COMM_WORLD);

  // Освобождение памяти
  delete [] pReceiveNum; 
  delete [] pReceiveInd;
}

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

Рассмотрим результаты вычислительных экспериментов, выполненных для оценки эффективности приведенного выше параллельного алгоритма умножения матрицы на вектор. Кроме того, используем полученные результаты для сравнения теоретических оценок и экспериментальных показателей времени вычислений и проверим тем самым точность полученных аналитических соотношений. Эксперименты проводились на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server (см. п. 1.2.3).

Определение параметров теоретических зависимостей (величин $$\tau , w, \alpha , \beta$$ ) осуществлялось следующим образом. Для оценки длительности $$\tau$$ базовой скалярной операции проводилось решение задачи умножения матрицы на вектор при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины $$\tau$$ было получено значение 1,93 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности $$\alpha$$ и пропускной способности $$\beta$$ соответственно 47 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.

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

Сравнение экспериментального времени $$T'_p$$ выполнения параллельного алгоритма и теоретического времени Tp, вычисленного в соответствии с выражением (6.8), представлено в таблице 6.2 и в графическом виде на рис. 6.3 и 6.4.

6.6. Умножение матрицы на вектор при разделении данных по столбцам

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

6.6.1. Определение подзадач и выделение информационных зависимостей

При таком способе разделения данных в качестве базовой подзадачи может быть выбрана операция умножения столбца матрицы А на один из элементов вектора b. Для организации вычислений в этом случае каждая базовая подзадача i, 0<=i<n, должна содержать i -й столбец матрицы А и i -е элементы bi и ci векторов b и с.

Параллельный алгоритм умножения матрицы на вектор начинается с того, что каждая базовая задача i выполняет умножение своего столбца матрицы А на элемент bi, в итоге в каждой подзадаче получается вектор c'(i) промежуточных результатов. Далее для получения элементов результирующего вектора с подзадачи должны обменяться своими промежуточными данными между собой (элемент j, 0<=j<n, частичного результата c'(i) подзадачи i, 0<=i<n, должен быть передан подзадаче j ). Данная обобщенная передача данных ( all-to-all communication или total exchange ) является наиболее общей коммуникационной процедурой и может быть реализована при помощи функции MPI_Alltoall библиотеки MPI. После выполнения передачи данных каждая базовая подзадача i, 0<=i<n, будет содержать n частичных значений c'i(j), 0<=j<n, сложением которых и определяется элемент ci вектора результата с (см. рис. 6.5).

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

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

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

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

Пусть, как и ранее, матрица А является квадратной, то есть m=n. На первом этапе вычислений каждый процессор умножает принадлежащие ему столбцы матрицы А на элементы вектора b, после умножения полученные значения суммируются для каждой строки матрицы А в отдельности$$c'_s(i)= \sum_{j=j_0}^{j_{l-1}} a_{sj} b_j, \; 0 \le s < n ,$$ ( j0 и jl-1 есть начальный и конечный индексы столбцов базовой подзадачи i, 0<=i<n). Поскольку размеры полосы матрицы А и блока вектора b равны n/p, то трудоемкость таких вычислений может оцениваться как T'= n2/p операций. После обмена данными между подзадачами на втором этапе вычислений каждый процессор суммирует полученные значения для своего блока результирующего вектора c. Количество суммируемых значений для каждого элемента ci вектора c совпадает с числом процессоров p, размер блока вектора c на одном процессоре равен n/p, и, тем самым, число выполняемых операций для второго этапа оказывается равным T''=n. С учетом полученных соотношений показатели ускорения и эффективности параллельного алгоритма могут быть выражены следующим образом:$$S_p = \frac{n^2}{n^2/p} = p, \quad E_p=\frac{n^2}{p \cdot (n^2/p)} = 1.$$

Теперь рассмотрим более точные соотношения для оценки времени выполнения параллельного алгоритма. С учетом ранее проведенных рассуждений время выполнения вычислительных операций алгоритма может быть оценено при помощи выражения$$T_p(cacl)=[n \cdot (2 \cdot \lceil n/p \rceil - 1) + n] \cdot \tau$$ (здесь, как и ранее, $$\tau$$ есть время выполнения одной элементарной скалярной операции).

Для выполнения операции обобщенной передачи данных рассмотрим два возможных способа реализации (см. также лекцию 3). Первый способ обеспечивается алгоритмом, согласно которому каждый процессор последовательно передает свои данные всем остальным процессорам вычислительном системы. Предположим, что процессоры могут одновременно отправлять и принимать сообщения и между любой парой процессоров имеется прямая линия связи, тогда оценка трудоемкости (время исполнения) такого алгоритма обобщенной передачи данных может быть определена как$$T_p^1(comm) = (p-1)(\alpha +w \lceil {n/p}\rceil / \beta).$$ (напомним, что $$\alpha$$ – латентность сети передачи данных, $$\beta$$ – пропускная способность сети, w – размер элемента данных в байтах).

Второй способ выполнения операции обмена данными рассмотрен в лекции 3, когда топология вычислительной сети может быть представлена в виде гиперкуба. Как было показано, выполнение такого алгоритма может быть осуществлено за $$\lceil log_2p\rceil$$ шагов, на каждом из которых каждый процессор передает и получает сообщение из n/2 элементов. Как результат, время операции передачи данных при таком подходе составляет величину:$$T_p^2(comm)=\lceil \log_2 p \rceil (\alpha + w(n/2)/ \beta)$$

С учетом (6.11) – (6.13) общее время выполнения параллельного алгоритма умножения матрицы на вектор при разбиении данных по столбцам выражается следующими соотношениями.

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

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

Сравнение экспериментального времени $$T_p^*$$ выполнения эксперимента и времени Tp, вычисленного по соотношениям (6.14), (6.15), представлено в таблице 6.4 и на рис. 6.6 и 6.7. Теоретическое время $$T_p^1$$ вычисляется согласно (6.14), а теоретическое время $$T_p^2$$ – в соответствии с (6.15).

6.7. Умножение матрицы на вектор при блочном разделении данных

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

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

Блочная схема разбиения матриц подробно рассмотрена в подразделе 6.1. При таком способе разделения данных исходная матрица A представляется в виде набора прямоугольных блоков:$$A = \begin{pmatrix} A_{00} A_{02} \ldots A_{0q-1} \\ \ldots \\ A_{s-11} A_{s-12} \ldots A_{s-1q-1} \end{pmatrix} ,$$ где Aij, 0<=i<s, 0<=j<q, есть блок матрицы:$$A = \begin{pmatrix} a_{i_0 j_0} a_{i_0 j_1} \ldots a_{i_0 j_{l-1}} \\ \ldots \\ a_{i_{k-1} j_0} a_{i_{k-1} j_1} \ldots a_{i_{k-1} j_{l-1}} \end{pmatrix} , \begin{aligned} i_{\nu} = ik+\nu, 0 \le \nu < k,\ k=m/s, \\ j_u = jl+u, 0 \le u \le l,\ l=n/q. \end{aligned}$$ (здесь, как и ранее, предполагается, что p=sxq, количество строк матрицы является кратным s, а количество столбцов – кратным q, то есть m=kxs и n=lxq ).

При использовании блочного представления матрицы A базовые подзадачи целесообразно определить на основе вычислений, выполняемых над матричными блоками. Для нумерации подзадач могут применяться индексы располагаемых в подзадачах блоков матрицы A, т.е. подзадача (i,j) содержит блок Aij. Помимо блока матрицы A каждая подзадача должна содержать также и блок вектора b. При этом для блоков одной и той же подзадачи должны соблюдаться определенные правила соответствия – операция умножения блока матрицы A ij может быть выполнена только, если блок вектора b'(i,j) имеет вид$$b'(i,j)=(b'_0(i,j),\ldots,b'_l-1(i,j)),\quad \text{где} \quad b'_u(i,j)=b_{j_u} , \; j_u=jl+u, \; 0\le u < l, \; l=n/q.$$

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

Рассмотрим общую схему параллельных вычислений для операции умножения матрицы на вектор при блочном разделении исходных данных. После перемножения блоков матрицы A и вектора b каждая подзадача (i,j) будет содержать вектор частичных результатов c'(i,j), определяемый в соответствии с выражениями$$c'_\nu(i,j)=\sum_{u=0}^{l-1} a_{i_{\nu} j_u} b_{j_u}, \; i_{\nu}=i+\nu, \; 0 \le \nu < k, \; k=m/s, \; j_u=jl+u, \; 0 \le u \le l, l=n/q$$

Поэлементное суммирование векторов частичных результатов для каждой горизонтальной полосы (данная процедура часто именуется операцией редукции – см. лекцию 3) блоков матрицы A позволяет получить результирующий вектор c$$c_{\eta} = \sum_{j=0}^{q-1} c'_{\nu}(i,j), \; 0 \le \eta < m, \; i=\eta/s, \; \nu = \eta - i \cdots .$$

Для размещения вектора c применим ту же схему, что и для исходного вектора b: организуем вычисления таким образом, чтобы при завершении расчетов вектор c располагался поблочно в каждой из вертикальных полос блоков матрицы A (тем самым, каждый блок вектора c должен быть продублирован по каждой горизонтальной полосе). Выполнение всех необходимых действий для этого – суммирование частичных результатов и дублирование блоков результирующего вектора — может быть обеспечено при помощи функции MPI_Allreduce библиотеки MPI.

Общая схема выполняемых вычислений для умножения матрицы на вектор при блочном разделении данных показана на рис. 6.8.

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

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

Размер блоков матрицы А может быть подобран таким образом, чтобы общее количество базовых подзадач совпадало с числом процессоров p. Так, например, если определить размер блочной решетки как p=sxq, то

k=m/s, l=n/q,

где k и l есть количество строк и столбцов в блоках матрицы А. Такой способ определения размера блоков приводит к тому, что объем вычислений в каждой подзадаче является равным, и, тем самым, достигается полная балансировка вычислительной нагрузки между процессорами.

Возможность выбора остается при определении размеров блочной структуры матрицы А. Большое количество блоков по горизонтали приводит к возрастанию числа итераций в операции редукции результатов блочного умножения, а увеличение размера блочной решетки по вертикали повышает объем передаваемых данных между процессорами. Простое, часто применяемое решение состоит в использовании одинакового количества блоков по вертикали и горизонтали, т.е.$$s=q=\sqrt{p}$$

Следует отметить, что блочная схема разделения данных является обобщением всех рассмотренных в данной лекции подходов. Действительно, при q=1 блоки сводятся к горизонтальным полосам матрицы А, при s=1 исходные данные разбиваются на вертикальные полосы.

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

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

Выполним анализ эффективности параллельного алгоритма умножения матрицы на вектор при обычных уже предположениях, что матрица А является квадратной, т.е. m=n. Будем предполагать также, что процессоры, составляющие многопроцессорную вычислительную систему, образуют прямоугольную решетку p=sxq ( s – количество строк в процессорной решетке, q – количество столбцов).

Общий анализ эффективности приводит к идеальным показателям параллельного алгоритма:$$S_p = \frac{n^2}{n^2/p} =p , \quad E_p = \frac{n^2}{p \cdot (n^2 / p)} = 1$$

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

Общее время умножения блоков матрицы А и вектора b может быть определено как$$T_p(calc) = \lceil n/s \rceil \cdot (2 \cdot \lceil n/q \rceil - 1) \cdot \tau .$$

Операция редукции данных может быть выполнена с использованием каскадной схемы и включает, тем самым, log2q итераций передачи сообщений размера $$w\lceil n/s \rceil$$. Как результат, оценка коммуникационных затрат параллельного алгоритма при использовании модели Хокни может быть определена при помощи следующего выражения$$T_p(comm)=(\alpha +w \lceil n/s \rceil / \beta) \lceil log_2 q \rceil$$

Таким образом, общее время выполнения параллельного алгоритма умножения матрицы на вектор при блочном разделении данных составляет$$T_p = \lceil n/s \rceil \cdot (2 \cdot \lceil n/q \rceil - 1) \cdot \tau + (\alpha + w \lceil n/s \rceil / \beta) \lceil \log_2 q \rceil .$$

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

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

Сравнение экспериментального времени $$T^*_p$$ выполнения эксперимента и теоретического времени T p, вычисленного в соответствии с выражением (6.19), представлено в таблице 6.5 и на рис. 6.10.

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

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

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

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

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

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

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

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

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

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