Анализ эффективности
Для анализа эффективности параллельных вычислений здесь и далее будут строиться два типа оценок. В первой из них трудоемкость алгоритмов оценивается в количестве вычислительных операций, необходимых для решения поставленной задачи, без учета затрат времени на передачу данных между процессорами, а длительность всех вычислительных операций считается одинаковой. Кроме того, константы в получаемых соотношениях, как правило, не указываются — для первого типа оценок важен прежде всего порядок сложности алгоритма, а не точное выражение времени выполнения вычислений. Как результат, в большинстве случаев подобные оценки получаются достаточно простыми и могут быть использованы для начального анализа эффективности разрабатываемых алгоритмов и методов.
Второй тип оценок направлен на формирование как можно более точных соотношений для предсказания времени выполнения алгоритмов. Получение таких оценок проводится, как правило, при помощи уточнения выражений, полученных на первом этапе. Для этого в имеющиеся соотношения вводятся параметры, задающие длительность выполнения операций, строятся оценки трудоемкости коммуникационных операций, указываются все необходимые константы. Точность получаемых выражений проверяется при помощи вычислительных экспериментов, по результатам которых время выполненных расчетов сравнивается с теоретически предсказанными оценками длительностей вычислений. Как результат, оценки подобного типа имеют, как правило, более сложный вид, но позволяют более точно оценивать эффективность разрабатываемых методов параллельных вычислений.
Рассмотрим трудоемкость алгоритма умножения матрицы на вектор. В случае если матрица А квадратная (m=n), последовательный алгоритм умножения матрицы на вектор имеет сложность T1=n2. В случае параллельных вычислений каждый процессор производит умножение только части (полосы) матрицы A на вектор b, размер этих полос равен n/p строк. При вычислении скалярного произведения одной строки матрицы и вектора необходимо произвести n операций умножения и (n-1) операций сложения.
Контрольные вопросы
Назовите основные способы распределения элементов матрицы между процессорами вычислительной системы.В чем состоит постановка задачи умножения матрицы на вектор?Какова вычислительная сложность последовательного алгоритма умножения матрицы на вектор?Почему при разработке параллельных алгоритмов умножения матрицы на вектор допустимо дублировать вектор-операнд на все процессоры?Какие подходы могут быть предложены для разработки параллельных алгоритмов умножения матрицы на вектор?Представьте общие схемы рассмотренных параллельных алгоритмов умножения матрицы на вектор.Проведите анализ и получите показатели эффективности для одного из рассмотренных алгоритмов.Какой из представленных алгоритмов умножения матрицы на вектор обладает лучшими показателями ускорения и эффективности? Может ли использование циклической схемы разделения данных повлиять на время работы каждого из представленных алгоритмов?Какие информационные взаимодействия выполняются для алгоритмов при ленточной схеме разделения данных? В чем различие необходимых операций передачи данных при разделении матрицы по строкам и столбцам?Какие информационные взаимодействия выполняются для блочного алгоритма умножения матрицы на вектор?Какая топология коммуникационной сети является целесообразной для каждого из рассмотренных алгоритмов?Дайте общую характеристику программной реализации алгоритма умножения матрицы на вектор при разделении данных по строкам. В чем могут состоять различия в программной реализации других рассмотренных алгоритмов?Какие функции библиотеки MPI оказались необходимыми при программной реализации алгоритмов?
Краткий обзор лекции
В данной лекции на примере задачи умножения матрицы на вектор рассматриваются возможные схемы разделения матриц между процессорами многопроцессорной вычислительной системы, которые могут быть использованы для организации параллельных вычислений при выполнении матричных операций. В числе излагаемых схем способы разбиения матриц на полосы (по вертикали или горизонтали) или на прямоугольные наборы элементов (блоки).
Далее в лекции с использованием рассмотренных способов разделения матриц подробно излагаются три возможных варианта параллельного выполнения операции умножения матрицы на вектор. Первый алгоритм основан на разделении матрицы между процессорами по строкам, второй – на разделении матрицы по столбцам, а третий – на блочном разделении данных. Каждый алгоритм представлен в соответствии с общей схемой, описанной в лекции 4, — вначале определяются базовые подзадачи, затем выделяются информационные зависимости подзадач, далее обсуждается масштабирование и распределение подзадач между процессорами. В завершение для каждого алгоритма проводится анализ эффективности получаемых параллельных вычислений и приводятся результаты вычислительных экспериментов. Для алгоритма умножения матрицы на вектор при ленточном разделении данных по строкам приводится возможный вариант программной реализации.
Полученные показатели эффективности показывают, что все используемые способы разделения данных приводят к равномерной балансировке вычислительной нагрузки, и отличия имеются только в трудоемкости выполняемых информационных взаимодействий между процессорами. В этом отношении представляется интересным проследить, как выбор способа разделения данных влияет на характер необходимых операций передачи данных, и выделить основные различия в коммуникационных действиях разных алгоритмов. Кроме того, важным является определение целесообразной структуры линий связи между процессорами для эффективного выполнения соответствующего параллельного алгоритма. Так, например, алгоритмы, основанные на ленточном разделении данных, ориентированы на топологию сети в виде гиперкуба или полного графа. Для реализации алгоритма, основанного на блочном разделении данных, необходимо наличие топологии решетки.
На рис. 6.11 на общем графике представлены показатели ускорения, полученные в результате выполнения вычислительных экспериментов для всех рассмотренных алгоритмов. Следует отметить, что дополнительные расчеты показывают , что при большем количестве процессоров и при большом размере матриц более эффективным становится блочный алгоритм умножения.
Рис. 6.11. Показатели ускорения рассмотренных параллельных алгоритмов умножения по результатам вычислительных экспериментов с матрицами размера 2000x2000 и векторами из 2000 элементов
Масштабирование и распределение подзадач по процессорам
В процессе умножения плотной матрицы на вектор количество вычислительных операций для получения скалярного произведения одинаково для всех базовых подзадач. Поэтому в случае когда число процессоров p меньше числа базовых подзадач m, мы можем объединить базовые подзадачи таким образом, чтобы каждый процессор выполнял несколько таких задач, соответствующих непрерывной последовательности строк матрицы А. В этом случае по окончании вычислений каждая базовая подзадача определяет набор элементов результирующего вектора с.
Распределение подзадач между процессорами вычислительной системы может быть выполнено произвольным образом.
Размер блоков матрицы А может быть подобран таким образом, чтобы общее количество базовых подзадач совпадало с числом процессоров p. Так, например, если определить размер блочной решетки как p=s·q, то
k=m/s, l=n/q,
где k и l есть количество строк и столбцов в блоках матрицы А. Такой способ определения размера блоков приводит к тому, что объем вычислений в каждой подзадаче является равным, и, тем самым, достигается полная балансировка вычислительной нагрузки между процессорами.
Возможность выбора остается при определении размеров блочной структуры матрицы А. Большое количество блоков по горизонтали приводит к возрастанию числа итераций в операции редукции результатов блочного умножения, а увеличение размера блочной решетки по вертикали повышает объем передаваемых данных между процессорами. Простое, часто применяемое решение состоит в использовании одинакового количества блоков по вертикали и горизонтали, т.е.
Следует отметить, что блочная схема разделения данных является обобщением всех рассмотренных в данной лекции подходов. Действительно, при q=1 блоки сводятся к горизонтальным полосам матрицы А, при s=1 исходные данные разбиваются на вертикальные полосы.
При решении вопроса распределения подзадач между процессорами должна учитываться возможность эффективного выполнения операции редукции. Возможный вариант подходящего способа распределения состоит в выделении для подзадач одной и той же горизонтальной полосы блоков множества процессоров, структура сети передачи данных между которыми имеет вид гиперкуба или полного графа.
Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью и равным объемом передаваемых данных. В случае когда количество столбцов матрицы превышает число процессоров, базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько соседних столбцов (в этом случае исходная матрица A разбивается на ряд вертикальных полос). При соблюдении равенства размера полос такой способ агрегации вычислений обеспечивает равномерность распределения вычислительной нагрузки по процессорам, составляющим многопроцессорную вычислительную систему.
Как и в предыдущем алгоритме, распределение подзадач между процессорами вычислительной системы может быть выполнено произвольным образом.
Обзор литературы
Задача умножения матрицы на вектор часто используется как демонстрационный пример параллельного программирования и, как результат, широко рассматривается в литературе. В качестве дополнительного учебного материала могут быть рекомендованы работы [[2], [51], [63]]. Широкое обсуждение вопросов параллельного выполнения матричных вычислений выполнено в работе [[30]].
При рассмотрении вопросов программной реализации параллельных методов может быть рекомендована работа [[23]]. В ней рассматривается хорошо известная и широко применяемая в практике параллельных вычислений программная библиотека численных методов ScaLAPACK.
Определение подзадач
Блочная схема разбиения матриц подробно рассмотрена в подразделе 6.1. При таком способе разделения данных исходная матрица A представляется в виде набора прямоугольных блоков:
где Aij, 0i<s, 0j<q, есть блок матрицы:
(здесь, как и ранее, предполагается, что p=s·q, количество строк матрицы является кратным s, а количество столбцов – кратным q, то есть m=k·s и n=l·q).
При использовании блочного представления матрицы A базовые подзадачи целесообразно определить на основе вычислений, выполняемых над матричными блоками. Для нумерации подзадач могут применяться индексы располагаемых в подзадачах блоков матрицы A, т.е. подзадача (i,j) содержит блок Aij. Помимо блока матрицы A каждая подзадача должна содержать также и блок вектора b. При этом для блоков одной и той же подзадачи должны соблюдаться определенные правила соответствия – операция умножения блока матрицы A ij может быть выполнена только, если блок вектора b'(i,j) имеет вид
Определение подзадач и выделение информационных зависимостей
При таком способе разделения данных в качестве базовой подзадачи может быть выбрана операция умножения столбца матрицы А на один из элементов вектора b. Для организации вычислений в этом случае каждая базовая подзадача i, 0i<n, должна содержать i-й столбец матрицы А и i-е элементы bi и ci векторов b и с.
Параллельный алгоритм умножения матрицы на вектор начинается с того, что каждая базовая задача i выполняет умножение своего столбца матрицы А на элемент bi, в итоге в каждой подзадаче получается вектор c'(i) промежуточных результатов. Далее для получения элементов результирующего вектора с подзадачи должны обменяться своими промежуточными данными между собой (элемент j, 0j<n, частичного результата
c'(i) подзадачи i, 0i<n, должен быть передан подзадаче j). Данная обобщенная передача данных (all-to-all communication или total exchange) является наиболее общей коммуникационной процедурой и может быть реализована при помощи функции MPI_Alltoall библиотеки MPI. После выполнения передачи данных каждая базовая подзадача i, 0i<n, будет содержать n частичных значений c'i(j), 0j<n, сложением которых и определяется элемент ci вектора результата с (см. рис. 6.5).
Рис. 6.5. Организация вычислений при выполнении параллельного алгоритма умножения матрицы на вектор с использованием разбиения матрицы по столбцам
Параллельные методы умножения матрицы на вектор
Матрицы и матричные операции широко используются при математическом моделировании самых разнообразных процессов, явлений и систем. Матричные вычисления составляют основу многих научных и инженерных расчетов – среди областей приложений могут быть указаны вычислительная математика, физика, экономика и др.
С учетом значимости эффективного выполнения матричных расчетов многие стандартные библиотеки программ содержат процедуры для различных матричных операций. Объем программного обеспечения для обработки матриц постоянно увеличивается – разрабатываются новые экономные структуры хранения для матриц специального типа (треугольных, ленточных, разреженных и т.п.), создаются различные высокоэффективные машинно-зависимые реализации алгоритмов, проводятся теоретические исследования для поиска более быстрых методов матричных вычислений.
Являясь вычислительно трудоемкими, матричные вычисления представляют собой классическую область применения параллельных вычислений. С одной стороны, использование высокопроизводительных многопроцессорных систем позволяет существенно повысить сложность решаемых задач. С другой стороны, в силу своей достаточно простой формулировки матричные операции предоставляют прекрасную возможность для демонстрации многих приемов и методов параллельного программирования.
В данной лекции обсуждаются методы параллельных вычислений для операции матрично-векторного умножения, в следующей лекции (лекция 7) излагается более общий случай – задача перемножения матриц. Важный вид матричных вычислений – решение систем линейных уравнений – представлен в лекции 8. Общий для всех перечисленных задач вопрос разделения обрабатываемых матриц между параллельно работающими процессорами рассматривается в первом подразделе лекции 6.
При изложении следующего материала будем полагать, что рассматриваемые матрицы являются плотными (dense), то есть число нулевых элементов в них незначительно по сравнению с общим количеством элементов матриц.
Последовательный алгоритм
Последовательный алгоритм умножения матрицы на вектор может быть представлен следующим образом.
Алгоритм 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).
Постановка задачи
В результате умножения матрицы А размерности m ? n и вектора b, состоящего из n элементов, получается вектор c размера m, каждый i-й элемент которого есть результат скалярного умножения i-й строки матрицы А (обозначим эту строчку ai) и вектора b:
|
(6.4) |
Тем самым получение результирующего вектора c предполагает повторение m однотипных операций по умножению строк матрицы A и вектора b. Каждая такая операция включает умножение элементов строки матрицы и вектора b (n операций) и последующее суммирование полученных произведений (n-1 операций). Общее количество необходимых скалярных операций есть величина
T1=m·(2n-1)
Принципы распараллеливания
Для многих методов матричных вычислений характерным является повторение одних и тех же вычислительных действий для разных элементов матриц. Данное свойство свидетельствует о наличии параллелизма по данным при выполнении матричных расчетов, и, как результат, распараллеливание матричных операций сводится в большинстве случаев к разделению обрабатываемых матриц между процессорами используемой вычислительной системы. Выбор способа разделения матриц приводит к определению конкретного метода параллельных вычислений; существование разных схем распределения данных порождает целый ряд параллельных алгоритмов матричных вычислений.
Наиболее общие и широко используемые способы разделения матриц состоят в разбиении данных на полосы (по вертикали или горизонтали) или на прямоугольные фрагменты (блоки).
1. Ленточное разбиение матрицы. При ленточном (block-striped) разбиении каждому процессору выделяется то или иное подмножество строк (rowwise или горизонтальное разбиение) или столбцов (columnwise или вертикальное разбиение) матрицы (рис. 6.1). Разделение строк и столбцов на полосы в большинстве случаев происходит на непрерывной (последовательной) основе. При таком подходе для горизонтального разбиения по строкам, например, матрица A представляется в виде (см. рис. 6.1)
|
(6.1) |
где ai=(ai1,ai2,...,ain), 0i<m, есть i-я строка матрицы A (предполагается, что количество строк m кратно числу процессоров p, т.е. m = k·p). Во всех алгоритмах матричного умножения и умножения матрицы на вектор, которые будут рассмотрены в этой и следующей лекциях, применяется разделение данных на непрерывной основе.
Другой возможный подход к формированию полос состоит в применении той или иной схемы чередования (цикличности) строк или столбцов. Как правило, для чередования используется число процессоров p – в этом случае при горизонтальном разбиении матрица A принимает вид
|
(6.2) |
Циклическая схема формирования полос может оказаться полезной для лучшей балансировки вычислительной нагрузки процессоров (например, при решении системы линейных уравнений с использованием метода Гаусса – см.
лекцию 8).
2. Блочное разбиение матрицы. При блочном ( chessboard block) разделении матрица делится на прямоугольные наборы элементов – при этом, как правило, используется разделение на непрерывной основе. Пусть количество процессоров составляет p = s·q, количество строк матрицы является кратным s, а количество столбцов – кратным q, то есть m = k·s и n = l·q. Представим исходную матрицу A в виде набора прямоугольных блоков следующим образом:
где Aij — блок матрицы, состоящий из элементов:
(6.3) |
При таком подходе целесообразно, чтобы вычислительная система имела физическую или, по крайней мере, логическую топологию процессорной решетки из s строк и q столбцов. В этом случае при разделении данных на непрерывной основе процессоры, соседние в структуре решетки, обрабатывают смежные блоки исходной матрицы. Следует отметить, однако, что и для блочной схемы может быть применено циклическое чередование строк и столбцов.
В данной лекции рассматриваются три параллельных алгоритма для умножения квадратной матрицы на вектор. Каждый подход основан на разном типе распределения исходных данных (элементов матрицы и вектора) между процессорами. Разделение данных меняет схему взаимодействия процессоров, поэтому каждый из представленных методов существенным образом отличается от двух остальных.
Рис. 6.1. Способы распределения элементов матрицы между процессорами вычислительной системы
Программная реализация
Представим возможный вариант параллельной программы умножения матрицы на вектор с использованием алгоритма разбиения матрицы по строкам. При этом реализация отдельных модулей не приводится, если их отсутствие не оказывает влияние на понимании общей схемы параллельных вычислений.
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; }
Разделение данных
При выполнении параллельных алгоритмов умножения матрицы на вектор, кроме матрицы А, необходимо разделить еще вектор b и вектор результата c. Элементы векторов можно продублировать, то есть скопировать все элементы вектора на все процессоры, составляющие многопроцессорную вычислительную систему, или разделить между процессорами. При блочном разбиении вектора из n элементов каждый процессор обрабатывает непрерывную последовательность из k элементов вектора (мы предполагаем, что размерность вектора n нацело делится на число процессоров, т.е. n= k·p).
Поясним, почему дублирование векторов b и
c между процессорами является допустимым решением (далее для простоты изложения будем полагать, что m=n). Векторы b и с состоят из n элементов, т.е. содержат столько же данных, сколько и одна строка или один столбец матрицы. Если процессор хранит строку или столбец матрицы и одиночные элементы векторов b и с, то общее число сохраняемых элементов имеет порядок O(n). Если процессор хранит строку (столбец) матрицы и все элементы векторов b и с, то общее число сохраняемых элементов также порядка O(n). Таким образом, при дублировании и при разделении векторов требования к объему памяти из одного класса сложности.
Результаты вычислительных экспериментов
Рассмотрим результаты вычислительных экспериментов, выполненных для оценки эффективности приведенного выше параллельного алгоритма умножения матрицы на вектор. Кроме того, используем полученные результаты для сравнения теоретических оценок и экспериментальных показателей времени вычислений и проверим тем самым точность полученных аналитических соотношений. Эксперименты проводились на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server (см. п. 1.2.3).
Определение параметров теоретических зависимостей (величин ?, w, , ?) осуществлялось следующим образом. Для оценки длительности ? базовой скалярной операции проводилось решение задачи умножения матрицы на вектор при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины ? было получено значение 1,93 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности и пропускной способности ? соответственно 47 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.
Результаты вычислительных экспериментов приведены в таблице 6.1. Эксперименты проводились с использованием двух, четырех и восьми процессоров. Времена выполнения алгоритмов указаны в секундах.
| 1000 | 0,0041 | 0,0021 | 1,8798 | 0,0017 | 2,4089 | 0,0175 | 0,2333 |
| 2000 | 0,016 | 0,0084 | 1,8843 | 0,0047 | 3,3388 | 0,0032 | 4,9443 |
| 3000 | 0,031 | 0,0185 | 1,6700 | 0,0097 | 3,1778 | 0,0059 | 5,1952 |
| 4000 | 0,062 | 0,0381 | 1,6263 | 0,0188 | 3,2838 | 0,0244 | 2,5329 |
| 5000 | 0,11 | 0,0574 | 1,9156 | 0,0314 | 3,4993 | 0,0150 | 7,3216 |
Результаты вычислительных экспериментов
Вычислительные эксперименты для оценки эффективности параллельного алгоритма умножения матрицы на вектор при разбиении данных по столбцам проводились при условиях, указанных в п. 6.5.5. Результаты вычислительных экспериментов приведены в таблице 6.3.
| 1000 | 0,0041 | 0,0022 | 1,8352 | 0,0015 | 3,1538 | 0,0008 | 4,9409 |
| 2000 | 0,016 | 0,0085 | 1,8799 | 0,0046 | 3,4246 | 0,0029 | 5,4682 |
| 3000 | 0,031 | 0,019 | 1,6315 | 0,0095 | 3,2413 | 0,0055 | 5,5456 |
| 4000 | 0,062 | 0,0331 | 1,8679 | 0,0168 | 3,6714 | 0,0090 | 6,8599 |
| 5000 | 0,11 | 0,0518 | 2,1228 | 0,0265 | 4,1361 | 0,0136 | 8,0580 |
Сравнение экспериментального времени
выполнения эксперимента и времени Tp, вычисленного по соотношениям (6.14), (6.15), представлено в таблице 6.4 и на рис. 6.6 и 6.7. Теоретическое время вычисляется согласно (6.14), а теоретическое время – в соответствии с (6.15).
Рис. 6.6. График зависимости теоретического и экспериментального времени выполнения параллельного алгоритма на четырех процессорах от объема исходных данных (ленточное разбиение матрицы по столбцам)
| 1000 | 0,0022 | 0,0021 | 0,0021 | 0,0013 | 0,0013 | 0,0014 | 0,0008 | 0,0011 | 0,0015 |
| 2000 | 0,0085 | 0,0082 | 0,0080 | 0,0046 | 0,0044 | 0,0044 | 0,0029 | 0,0027 | 0,0031 |
| 3000 | 0,019 | 0,0177 | 0,0177 | 0,0095 | 0,0094 | 0,0094 | 0,0055 | 0,0054 | 0,0056 |
| 4000 | 0,0331 | 0,0313 | 0,0313 | 0,0168 | 0,0163 | 0,0162 | 0,0090 | 0,0090 | 0,0091 |
| 5000 | 0,0518 | 0,0487 | 0,0487 | 0,0265 | 0,0251 | 0,0251 | 0,0136 | 0,0135 | 0,0136 |
Рис. 6.7. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма умножения матрицы на вектор (ленточное разбиение матрицы по столбцам) для разных размеров матриц
Результаты вычислительных экспериментов
Вычислительные эксперименты для оценки эффективности параллельного алгоритма проводились при тех же условиях, что и ранее выполненные расчеты (см. п. 6.5.5). Результаты экспериментов приведены в таблице 6.5. Вычисления проводились с использованием четырех и девяти процессоров.
Сравнение экспериментального времени
выполнения эксперимента и теоретического времени T
p, вычисленного в соответствии с выражением (6.19), представлено в таблице 6.5 и на рис. 6.10.
| 1000 | 0,0041 | 0,0028 | 1,4260 | 0,0011 | 3,7998 |
| 2000 | 0,016 | 0,0099 | 1,6127 | 0,0095 | 3,2614 |
| 3000 | 0,031 | 0,0214 | 1,4441 | 0,0095 | 3,2614 |
| 4000 | 0,062 | 0,0381 | 1,6254 | 0,0175 | 3,5420 |
| 5000 | 0,11 | 0,0583 | 1,8860 | 0,0263 | 4,1755 |
Рис. 6.9. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма умножения матрицы на вектор (блочное разбиение матрицы) для разных размеров матриц
| 1000 | 0,0025 | 0,0028 | 0,0012 | 0,0011 |
| 2000 | 0,0095 | 0,0099 | 0,0043 | 0,0042 |
| 3000 | 0,0212 | 0,0214 | 0,0095 | 0,0095 |
| 4000 | 0,0376 | 0,0381 | 0,0168 | 0,0175 |
| 5000 | 0,0586 | 0,0583 | 0,0262 | 0,0263 |
Рис. 6.10. График зависимости экспериментального и теоретического времени проведения эксперимента на четырех процессорах от объема исходных данных (блочное разбиение матрицы)
Умножение матрицы на вектор при блочном разделении данных
Рассмотрим теперь параллельный алгоритм умножения матрицы на вектор, который основан на ином способе разделения данных – на разбиении матрицы на прямоугольные фрагменты (блоки).
Умножение матрицы на вектор при разделении данных по строкам
Рассмотрим в качестве первого примера организации параллельных матричных вычислений алгоритм умножения матрицы на вектор, основанный на представлении матрицы непрерывными наборами (горизонтальными полосами) строк. При таком способе разделения данных в качестве базовой подзадачи может быть выбрана операция скалярного умножения одной строки матрицы на вектор.
Рассмотрим теперь другой подход к параллельному умножению матрицы на вектор, основанный на разделении исходной матрицы на непрерывные наборы (вертикальные полосы) столбцов.
Выделение информационных зависимостей
Рассмотрим общую схему параллельных вычислений для операции умножения матрицы на вектор при блочном разделении исходных данных. После перемножения блоков матрицы A и вектора b каждая подзадача (i,j) будет содержать вектор частичных результатов c'(i,j), определяемый в соответствии с выражениями
Поэлементное суммирование векторов частичных результатов для каждой горизонтальной полосы (данная процедура часто именуется операцией редукции – см. лекцию 3) блоков матрицы A позволяет получить результирующий вектор c
Для размещения вектора c применим ту же схему, что и для исходного вектора b: организуем вычисления таким образом, чтобы при завершении расчетов вектор c располагался поблочно в каждой из вертикальных полос блоков матрицы A (тем самым, каждый блок вектора c должен быть продублирован по каждой горизонтальной полосе). Выполнение всех необходимых действий для этого – суммирование частичных результатов и дублирование блоков результирующего вектора — может быть обеспечено при помощи функции MPI_Allreduce библиотеки MPI.
Общая схема выполняемых вычислений для умножения матрицы на вектор при блочном разделении данных показана на рис. 6.8.
Рис. 6.8. Общая схема параллельного алгоритма умножения матрицы на вектор при блочном разделении данных: a) исходное распределение результатов, б) распределение векторов частичных результатов, в) распределение блоков результирующего вектора c
Рассмотрев представленную схему параллельных вычислений, можно сделать вывод, что информационная зависимость базовых подзадач проявляется только на этапе суммирования результатов перемножения блоков матрицы A и блоков вектора b. Выполнение таких расчетов может быть выполнено по обычной каскадной схеме (см. лекцию 3), и, как результат, характер имеющихся информационных связей для подзадач одной и той же горизонтальной полосы блоков соответствует структуре двоичного дерева.
Для выполнения базовой подзадачи скалярного произведения процессор должен содержать соответствующую строку матрицы А и копию вектора b. После завершения вычислений каждая базовая подзадача определяет один из элементов вектора результата c.
Для объединения результатов расчета и получения полного вектора c на каждом из процессоров вычислительной системы необходимо выполнить операцию обобщенного сбора данных (см. лекцию 4), в которой каждый процессор передает свой вычисленный элемент вектора c всем остальным процессорам. Этот шаг можно выполнить, например, с использованием функции MPI_Allgather из библиотеки MPI.
В общем виде схема информационного взаимодействия подзадач в ходе выполняемых вычислений показана на рис. 6.2.
Рис. 6.2. Организация вычислений при выполнении параллельного алгоритма умножения матрицы на вектор, основанного на разделении матрицы по строкам
Выполните реализацию параллельного алгоритма, основанного
Выполните реализацию параллельного алгоритма, основанного на ленточном разбиении матрицы на вертикальные полосы. Постройте ие оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты реальных экспериментов с ранее подготовленными теоретическими оценками.Выполните реализацию параллельного алгоритма, основанного на разбиении матрицы на блоки. Постройте теоретические оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты реальных экспериментов с ранее подготовленными теоретическими оценками.
Алгоритм Фокса умножения матриц при блочном разделении данных
При построении параллельных способов выполнения матричного умножения наряду с рассмотрением матриц в виде наборов строк и столбцов широко используется блочное представление матриц. Рассмотрим более подробно данный способ организации вычислений.
Алгоритм Кэннона умножения матриц при блочном разделении данных
Рассмотрим еще один параллельный алгоритм матричного умножения, основанный на блочном разбиении матриц, – алгоритм Кэннона (Cannon).
Анализ эффективности
Выполним анализ эффективности первого параллельного алгоритма умножения матриц.
Общая трудоемкость последовательного алгоритма, как уже отмечалось ранее, является пропорциональной n3. Для параллельного алгоритма на каждой итерации каждый процессор выполняет умножение имеющихся на процессоре полос матрицы А и матрицы В (размер полос равен n/p, и, как результат, общее количество выполняемых при этом умножении операций равно n3/p2). Поскольку число итераций алгоритма совпадает с количеством процессоров, сложность параллельного алгоритма без учета затрат на передачу данных может быть определена при помощи выражения
|
(7.4) |
С учетом этой оценки показатели ускорения и эффективности данного параллельного алгоритма матричного умножения принимают вид:
|
(7.5) |
Таким образом, общий анализ сложности дает идеальные показатели эффективности параллельных вычислений. Для уточнения полученных соотношений оценим более точно количество вычислительных операций алгоритма и учтем затраты на выполнение операций передачи данных между процессорами.
С учетом числа и длительности выполняемых операций время выполнения вычислений параллельного алгоритма может быть оценено следующим образом:
|
(7.6) |
(здесь, как и ранее, ? есть время выполнения одной элементарной скалярной операции).
Для оценки коммуникационной сложности параллельных вычислений будем предполагать, что все операции передачи данных между процессорами в ходе одной итерации алгоритма могут быть выполнены параллельно. Объем передаваемых данных между процессорами определяется размером полос и составляет n/p строк или столбцов длины n. Общее количество параллельных операций передачи сообщений на единицу меньше числа итераций алгоритма (на последней итерации передача данных не является обязательной). Тем самым, оценка трудоемкости выполняемых операций передачи данных может быть определена как
|
(7.7) |
где – латентность, ? – пропускная способность сети передачи данных, а w есть размер элемента матрицы в байтах.
С учетом полученных соотношений общее время выполнения параллельного алгоритма матричного умножения определяется следующим выражением:
|
(7.8) |
Определим вычислительную сложность данного алгоритма Фокса. Построение оценок будет происходить при условии выполнения всех ранее выдвинутых предположений: все матрицы являются квадратными размера n?n, количество блоков по горизонтали и вертикали являются одинаковым и равным q (т.е. размер всех блоков равен k?k, k=n/q), процессоры образуют квадратную решетку и их количество равно p=q2.
Как уже отмечалось, алгоритм Фокса требует для своего выполнения q итераций, в ходе которых каждый процессор перемножает свои текущие блоки матриц А и В и прибавляет результаты умножения к текущему значению блока матрицы C. С учетом выдвинутых предположений общее количество выполняемых при этом операций будет иметь порядок n3/p. Как результат, показатели ускорения и эффективности алгоритма имеют вид:
|
(7.9) |
Общий анализ сложности снова дает идеальные показатели эффективности параллельных вычислений. Уточним полученные соотношения — для этого укажем более точно количество вычислительных операций алгоритма и учтем затраты на выполнение операций передачи данных между процессорами.
Определим количество вычислительных операций. Сложность выполнения скалярного умножения строки блока матрицы A на столбец блока матрицы В можно оценить как 2(n/q)-1. Количество строк и столбцов в блоках равно n/q и, как результат, трудоемкость операции блочного умножения оказывается равной (n2/p)(2n/q-1). Для сложения блоков требуется n2/p операций. С учетом всех перечисленных выражений время выполнения вычислительных операций алгоритма Фокса может быть оценено следующим образом:
|
(7.10) |
(напомним, что ? есть время выполнения одной элементарной скалярной операции).
Оценим затраты на выполнение операций передачи данных между процессорами. На каждой итерации алгоритма перед умножением блоков один из процессоров строки процессорной решетки рассылает свой блок матрицы A остальным процессорам своей строки. Как уже отмечалось ранее, при топологии сети в виде гиперкуба или полного графа выполнение этой операции может быть обеспечено за log2q шагов, а объем передаваемых блоков равен n2/p.
Перед проведением анализа эффективности следует отметить, что алгоритм Кэннона отличается от метода Фокса только видом выполняемых в ходе вычислений коммуникационных операций. Как результат, используя оценки времени выполнения вычислительных операций, приведенные в п. 7.4.4, проведем только анализ коммуникационной сложности алгоритма Кэннона.
В соответствии с правилами алгоритма на этапе инициализации производится перераспределение блоков матриц А и B при помощи циклического сдвига матричных блоков по строкам и столбцам процессорной решетки. Трудоемкость выполнения такой операции передачи данных существенным образом зависит от топологии сети. Для сети со структурой полного графа все необходимые пересылки блоков могут быть выполнены одновременно (т.е. длительность операции оказывается равной времени передачи одного матричного блока между соседними процессорами). Для сети с топологией гиперкуба операция циклического сдвига может потребовать выполнения log2q итераций. Для сети с кольцевой структурой связей необходимое количество итераций оказывается равным q-1 – более подробно методы выполнения операции циклического сдвига рассмотрены в лекции 3. Используем для построения оценки коммуникационной сложности этапа инициализации вариант топологии полного графа как более соответствующего кластерным вычислительным системам, время выполнения начального перераспределения блоков может оцениваться как
|
(7.14) |
(выражение n2/p определяет размер пересылаемых блоков, а коэффициент 2 соответствует двум выполняемым операциям циклического сдвига).
Оценим теперь затраты на передачу данных между процессорами при выполнении основной части алгоритма Кэннона. На каждой итерации алгоритма после умножения матричных блоков процессоры передают свои блоки предыдущим процессорам по строкам (для блоков матрицы A) и столбцам (для блоков матрицы B) процессорной решетки. Эти операции также могут быть выполнены процессорами параллельно, и, тем самым, длительность таких коммуникационных действий составляет:
|
() |
Поскольку количество итераций алгоритма Кэннона равно q, то с учетом оценки (7.10) общее время выполнения параллельных вычислений может быть определено при помощи следующего соотношения:
|
(7.16) |
(в используемых выражениях параметр определяет размер процессорной решетки).
Контрольные вопросы
В чем состоит постановка задачи умножения матриц? Приведите примеры задач, в которых используется операция умножения матриц.Приведите примеры различных последовательных алгоритмов выполнения операции умножения матриц. Отличается ли их вычислительная трудоемкость?Какие способы разделения данных используются при разработке параллельных алгоритмов матричного умножения?Представьте общие схемы рассмотренных параллельных алгоритмов умножения матриц.Проведите анализ и получите показатели эффективности ленточного алгоритма при горизонтальном разбиении перемножаемых матриц.Какие информационные взаимодействия выполняются для алгоритмов при ленточной схеме разделения данных? Какие информационные взаимодействия выполняются для блочных алгоритмов умножения матриц?Какая топология коммуникационной сети является целесообразной для каждого из рассмотренных алгоритмов?Какой из рассмотренных алгоритмов характеризуется наименьшими и наибольшими требованиями к объему необходимой памяти?Какой из рассмотренных алгоритмов обладает наилучшими показателями ускорения и эффективности?Оцените возможность выполнения матричного умножения как последовательности операций умножения матрицы на вектор.Дайте общую характеристику программной реализации алгоритма Фокса. В чем могут состоять различия в программной реализации других рассмотренных алгоритмов?Какие функции библиотеки MPI оказались необходимыми при программной реализации алгоритмов?
Краткий обзор лекции
В данной лекции рассмотрены три параллельных метода для выполнения операции матричного умножения. Первый алгоритм основан на ленточном разделении матриц между процессорами. В лекции приведены два различных варианта этого алгоритма. Первый вариант алгоритма основан на различном разделении перемножаемых матриц – первая матрица (матрица A) разбивается на горизонтальные полосы, а вторая матрица (матрица B) делится на вертикальные полосы. Второй вариант ленточного алгоритма использует разбиение обеих матриц на горизонтальные полосы.
Далее в лекции рассматриваются широко известные алгоритмы Фокса и Кэннона, основанные на блочном разделении матриц. При использовании одинаковой схемы разбиения матриц данные алгоритмы отличаются характером выполняемых операций передачи данных. Для алгоритма Фокса в ходе вычислений осуществляется рассылка и циклический сдвиг блоков матриц, в алгоритме Кэннона выполняется только операция циклического сдвига.
Различие в способах разбиения данных приводит к разным топологиям коммуникационной сети, при которых выполнение параллельных алгоритмов является наиболее эффективным. Так, алгоритмы, основанные на ленточном разделении данных, ориентированы на топологию сети в виде гиперкуба или полного графа. Для реализации алгоритмов, основанных на блочном разделении данных, необходимо наличие топологии решетки.
На рис. 7.12 на общем графике представлены показатели ускорения, полученные в результате выполнения вычислительных экспериментов для всех рассмотренных алгоритмов. Выполненные расчеты показывают, что при большем количестве процессоров более эффективными становятся блочные алгоритмы умножения матриц.
Рис. 7.12. Ускорение трех параллельных алгоритмов при умножении матриц с использованием четырех процессоров
Масштабирование и распределение подзадач по процессорам
Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью и равным объемом передаваемых данных. Когда размер матриц n оказывается больше, чем число процессоров p, базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько соседних строк и столбцов перемножаемых матриц. В этом случае исходная матрица A разбивается на ряд горизонтальных полос, а матрица B представляется в виде набора вертикальных (для первого алгоритма) или горизонтальных (для второго алгоритма) полос. Размер полос при этом следует выбрать равным k=n/p (в предположении, что n кратно p), что позволит по-прежнему обеспечить равномерность распределения вычислительной нагрузки по процессорам, составляющим многопроцессорную вычислительную систему.
Для распределения подзадач между процессорами может быть использован любой способ, обеспечивающий эффективное представление кольцевой структуры информационного взаимодействия подзадач. Для этого достаточно, например, чтобы подзадачи, являющиеся соседними в кольцевой топологии, располагались на процессорах, между которыми имеются прямые линии передачи данных.
В рассмотренной схеме параллельных вычислений количество блоков может варьироваться в зависимости от выбора размера блоков – эти размеры могут быть подобраны таким образом, чтобы общее количество базовых подзадач совпадало с числом процессоров p. Так, например, в наиболее простом случае, когда число процессоров представимо в виде p=2 (т.е. является полным квадратом), можно выбрать количество блоков в матрицах по вертикали и горизонтали равным (т.е. q=). Такой способ определения количества блоков приводит к тому, что объем вычислений в каждой подзадаче является одинаковым и тем самым достигается полная балансировка вычислительной нагрузки между процессорами. В более общем случае при произвольных количестве процессоров и размерах матриц балансировка вычислений может отличаться от абсолютно одинаковой, но, тем не менее, при надлежащем выборе параметров может быть распределена между процессорами равномерно в рамках требуемой точности.
Для эффективного выполнения алгоритма Фокса, в котором базовые подзадачи представлены в виде квадратной решетки и в ходе вычислений выполняются операции передачи блоков по строкам и столбцам решетки подзадач, наиболее адекватным решением является организация множества имеющихся процессоров также в виде квадратной решетки. В этом случае можно осуществить непосредственное отображение набора подзадач на множество процессоров – базовую подзадачу (i,j) следует располагать на процессоре Pi,j. Необходимая структура сети передачи данных может быть обеспечена на физическом уровне, если топология вычислительной системы имеет вид решетки или полного графа.
Как и ранее в методе Фокса, для алгоритма Кэннона размер блоков может быть подобран таким образом, чтобы количество базовых подзадач совпадало с числом имеющихся процессоров. Поскольку объемы вычислений в каждой подзадаче являются равными, это обеспечивает полную балансировку вычислительной нагрузки между процессорами.
Для распределения подзадач между процессорами может быть применен подход, использованный в алгоритме Фокса, – множество имеющихся процессоров представляется в виде квадратной решетки и размещение базовых подзадач (i,j) осуществляется на процессорах Pi,j соответствующих узлов процессорной решетки. Необходимая структура сети передачи данных, как и ранее, может быть обеспечена на физическом уровне при топологии вычислительной системы в виде решетки или полного графа.
Обзор литературы
Задача умножения матриц широко рассматривается в литературе. В качестве дополнительного учебного материала могут быть рекомендованы работы [[2], [51], [63]]. Широкое обсуждение вопросов параллельного выполнения матричных вычислений приводится в работе [[30]].
При рассмотрении вопросов программной реализации параллельных методов может быть рекомендована работа [[23]]. В ней рассматривается хорошо известная и широко используемая в практике параллельных вычислений программная библиотека численных методов ScaLAPACK.
Определение подзадач
Из определения операции матричного умножения следует, что вычисление всех элементов матрицы С может быть выполнено независимо друг от друга. Как результат, возможный подход для организации параллельных вычислений состоит в использовании в качестве базовой подзадачи процедуры определения одного элемента результирующей матрицы С. Для проведения всех необходимых вычислений каждая подзадача должна содержать по одной строке матрицы А и одному столбцу матрицы В. Общее количество получаемых при таком подходе подзадач оказывается равным n2
(по числу элементов матрицы С).
Рассмотрев предложенный подход, можно отметить, что достигнутый уровень параллелизма является в большинстве случаев избыточным. Обычно при проведении практических расчетов такое количество сформированных подзадач превышает число имеющихся процессоров и делает неизбежным этап укрупнения базовых задач. В этом плане может оказаться полезной агрегация вычислений уже на шаге выделения базовых подзадач. Возможное решение может состоять в объединении в рамках одной подзадачи всех вычислений, связанных не с одним, а с несколькими элементами результирующей матрицы С. Для дальнейшего рассмотрения определим базовую задачу как процедуру вычисления всех элементов одной из строк матрицы С. Такой подход приводит к снижению общего количества подзадач до величины n.
Для выполнения всех необходимых вычислений базовой подзадаче должны быть доступны одна из строк матрицы A и все столбцы матрицы B. Простое решение этой проблемы – дублирование матрицы B во всех подзадачах – является, как правило, неприемлемым в силу больших затрат памяти для хранения данных. Поэтому организация вычислений должна быть построена таким образом, чтобы в каждый текущий момент времени подзадачи содержали лишь часть данных, необходимых для проведения расчетов, а доступ к остальной части данных обеспечивался бы при помощи передачи данных между процессорами. Два возможных способа выполнения параллельных вычислений подобного типа рассмотрены далее в п. 7.3.2.
Как и при рассмотрении алгоритма Фокса, в качестве базовой подзадачи выберем вычисления, связанные с определением одного из блоков результирующей матрицы C. Как уже отмечалось ранее, для вычисления элементов этого блока подзадача должна иметь доступ к элементам горизонтальной полосы матрицы А и элементам вертикальной полосы матрицы В.
Блочная схема разбиения матриц подробно изложена в первом разделе лекции 6. При таком способе разделения данных исходные матрицы А, В и результирующая матрица С представляются в виде наборов блоков. Для более простого изложения следующего материала будем предполагать далее, что все матрицы являются квадратными размера n?n, количество блоков по горизонтали и вертикали одинаково и равно q (т.е. размер всех блоков равен k?k, k=n/q). При таком представлении данных операция матричного умножения матриц А и B в блочном виде может быть представлена так:
где каждый блок Cij матрицы C определяется в соответствии с выражением
При блочном разбиении данных для определения базовых подзадач естественным представляется взять за основу вычисления, выполняемые над матричными блоками. С учетом сказанного определим базовую подзадачу как процедуру вычисления всех элементов одного из блоков матрицы С.
Для выполнения всех необходимых вычислений базовым подзадачам должны быть доступны соответствующие наборы строк матрицы A и столбцов матрицы B. Размещение всех требуемых данных в каждой подзадаче неизбежно приведет к дублированию и к значительному росту объема используемой памяти. Как результат, вычисления должны быть организованы таким образом, чтобы в каждый текущий момент времени подзадачи содержали лишь часть необходимых для проведения расчетов данных, а доступ к остальной части данных обеспечивался бы при помощи передачи данных между процессорами. Один из возможных подходов – алгоритм Фокса (Fox) – рассмотрен далее в данном подразделе. Второй способ – алгоритм Кэннона (Cannon) – приводится в подразделе 7.5.
Последовательный алгоритм
Последовательный алгоритм умножения матриц представляется тремя вложенными циклами:
Алгоритм 7.1. Последовательный алгоритм умножения двух квадратных матриц
// Алгоритм 7.1 // Последовательный алгоритм умножения матриц double MatrixA[Size][Size]; double MatrixB[Size][Size]; double MatrixC[Size][Size]; int i,j,k; ... for (i=0; i<Size; i++){ for (j=0; j<Size; j++){ MatrixC[i][j] = 0; for (k=0; k<Size; k++){ MatrixC[i][j] = MatrixC[i][j] + MatrixA[i][k]*MatrixB[k][j]; } } }
Этот алгоритм является итеративным и ориентирован на последовательное вычисление строк матрицы С. Действительно, при выполнении одной итерации внешнего цикла (цикла по переменной i) вычисляется одна строка результирующей матрицы (см. рис. 7.1).
Рис. 7.1. На первой итерации цикла по переменной i используется первая строка матрицы A и все столбцы матрицы B для того, чтобы вычислить элементы первой строки результирующей матрицы С
Поскольку каждый элемент результирующей матрицы есть скалярное произведение строки и столбца исходных матриц, то для вычисления всех элементов матрицы С размером n?n необходимо выполнить n2·(2n–1) скалярных операций и затратить время
|
(7.3) |
где ? есть время выполнения одной элементарной скалярной операции.
Постановка задачи
Умножение матрицы A размера m?n и матрицы B размера n?l приводит к получению матрицы С размера m?l, каждый элемент которой определяется в соответствии с выражением:
|
(7.1) |
Как следует из (7.1), каждый элемент результирующей матрицы С есть скалярное произведение соответствующих строки матрицы A и столбца матрицы B:
|
(7.2) |
Этот алгоритм предполагает выполнение m·n·l операций умножения и столько же операций сложения элементов исходных матриц. При умножении квадратных матриц размера n?n количество выполненных операций имеет порядок O(n3). Известны последовательные алгоритмы умножения матриц, обладающие меньшей вычислительной сложностью (например, алгоритм Страссена (Strassen’s algorithm)), но эти алгоритмы требуют больших усилий для их освоения, и поэтому в данной лекции при разработке параллельных методов в качестве основы будет использоваться приведенный выше последовательный алгоритм. Также будем предполагать далее, что все матрицы являются квадратными и имеют размер n?n.
Условия выполнения программы: все матрицы
|
// Программа 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(); } |
| Пример 7.1. |
| Закрыть окно |
и инициализации исходных данных void
|
// Функция для выделения памяти и инициализации исходных данных 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); } } |
| Пример 7.3. |
| Закрыть окно |
// Функция для выделения памяти и инициализации исходных данных
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
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);
}
}
Программная реализация
Представим возможный вариант программной реализации алгоритма Фокса для умножения матриц при блочном представлении данных. Приводимый программный код содержит основные модули параллельной программы, отсутствие отдельных вспомогательных функций не сказывается на общем понимании реализуемой схемы параллельных вычислений.
1. Главная функция программы. Определяет основную логику работы алгоритма, последовательно вызывает необходимые подпрограммы.
Пример 7.1.
(html, txt)
2. Функция CreateGridCommunicators. Данная функция создает коммуникатор в виде двумерной квадратной решетки, определяет координаты каждого процесса в этой решетке, а также создает коммуникаторы отдельно для каждой строки и каждого столбца.
Создание решетки производится при помощи функции MPI_Cart_create (вектор Periodic определяет возможность передачи сообщений между граничными процессами строк и столбцов создаваемой решетки). После создания решетки каждый процесс параллельной программы будет иметь координаты своего положения в решетке; получение этих координат обеспечивается при помощи функции MPI_Cart_coords.
Формирование топологий завершается созданием множества коммуникаторов для каждой строки и каждого столбца решетки в отдельности (функция MPI_Cart_sub).
Пример 7.2.
(html, txt)
3. Функция ProcessInitialization. Данная функция определяет параметры решаемой задачи (размеры матриц и их блоков), выделяет память для хранения данных и осуществляет ввод исходных матриц (или формирует их при помощи какого-либо датчика случайных чисел). Всего в каждом процессе должна быть выделена память для хранения четырех блоков – для указателей на выделенную память используются переменные pAblock, pBblock, pCblock, pMatrixAblock. Первые три указателя определяют блоки матриц A, B и C соответственно. Следует отметить, что содержимое блоков pAblock и pBblock постоянно меняется в соответствии с пересылкой данных между процессами, в то время как блок pMatrixAblock матрицы A остается неизменным и применяется при рассылках блоков по строкам решетки процессов (см.
// Циклический сдвиг блоков матрицы В вдоль столбца процессной // решетки 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); } }
Результаты вычислительных экспериментов
Эксперименты проводились на вычислительном кластере на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server.
Для оценки длительности ? базовой скалярной операции проводилось решение задачи умножения матриц при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины ? было получено значение 6,4 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности a и пропускной способности b соответственно 130 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.
Результаты вычислительных экспериментов приведены в таблице 7.1. Эксперименты выполнялись с использованием двух, четырех и восьми процессоров.
| 500 | 0,8752 | 0,3758 | 2,3287 | 0,1535 | 5,6982 | 0,0968 | 9,0371 |
| 1000 | 12,8787 | 5,4427 | 2,3662 | 2,2628 | 5,6912 | 0,6998 | 18,4014 |
| 1500 | 43,4731 | 20,9503 | 2,0750 | 11,0804 | 3,9234 | 5,1766 | 8,3978 |
| 2000 | 103,0561 | 45,7436 | 2,2529 | 21,6001 | 4,7710 | 9,4127 | 10,9485 |
| 2500 | 201,2915 | 99,5097 | 2,0228 | 56,9203 | 3,5363 | 18,3303 | 10,9813 |
| 3000 | 347,8434 | 171,9232 | 2,0232 | 111,9642 | 3,1067 | 45,5482 | 7,6368 |
Рис. 7.4. Зависимость ускорения от количества процессоров при выполнении первого параллельного алгоритма матричного умножения при ленточной схеме распределения данных
Сравнение экспериментального времени выполнения эксперимента и теоретического времени Tp из формулы (7.8) представлено в таблице 7.2 и на рис. 7.5.
Рис. 7.5. График зависимости от объемап исходных данных теоретического и экспериментального времени выполнения параллельного алгоритма на двух процессорах (ленточная схема разбиения данных
Вычислительные эксперименты для оценки эффективности параллельного алгоритма проводились при тех же условиях, что и ранее выполненные (см. п. 7.3.5). Результаты экспериментов с использованием четырех и девяти процессоров приведены в таблице 7.3.
| 500 | 0,8527 | 0,2190 | 3,8925 | 0,1468 | 5,8079 |
| 1000 | 12,8787 | 3,0910 | 4,1664 | 2,1565 | 5,9719 |
| 1500 | 43,4731 | 10,8678 | 4,0001 | 7,2502 | 5,9960 |
| 2000 | 103,0561 | 24,1421 | 4,2687 | 21,4157 | 4,8121 |
| 2500 | 201,2915 | 51,4735 | 3,9105 | 41,2159 | 4,8838 |
| 3000 | 347,8434 | 87,0538 | 3,9957 | 58,2022 | 5,9764 |
Рис. 7.7. Зависимость ускорения от размера матриц при выполнении параллельного алгоритма Фокса
Рис. 7.8. График зависимости экспериментального и теоретического времени выполнения алгоритма Фокса на четырех процессорах
| 500 | 0,4217 | 0,2190 | 0,2200 | 0,1468 |
| 1000 | 3,2970 | 3,0910 | 1,5924 | 2,1565 |
| 1500 | 11,0419 | 10,8678 | 5,1920 | 7,2502 |
| 2000 | 26,0726 | 24,1421 | 12,0927 | 21,4157 |
| 2500 | 50,8049 | 51,4735 | 23,3682 | 41,2159 |
| 3000 | 87,6548 | 87,0538 | 40,0923 | 58,2022 |
Сравнение времени выполнения эксперимента и теоретического времени Tp, вычисленного в соответствии с выражением (7.13), представлено в таблице 7.4 и на рис. 7.8.
Вычислительные эксперименты для оценки эффективности параллельного алгоритма проводились при тех же условиях, что и ранее выполненные (см. п. 7.3.5). Результаты экспериментов для случаев четырех и девяти процессоров приведены в таблице 7.5.
| 1000 | 12,8787 | 3,0806 | 4,1805 | 1,1889 | 10,8324 |
| 1500 | 43,4731 | 11,1716 | 3,8913 | 4,6310 | 9,3872 |
| 2000 | 103,0561 | 24,0502 | 4,2850 | 14,4759 | 7,1191 |
| 2500 | 201,2915 | 53,1444 | 3,7876 | 23,5398 | 8,5511 |
| 3000 | 347,8434 | 88,2979 | 3,9394 | 36,3688 | 9,5643 |
Сравнение времени выполнения эксперимента и теоретического времени Tp, вычисленного в соответствии с выражением (7.16), представлено в таблице 7.6 и на рис. 7.11.
| 1000 | 3,4485 | 3,0806 | 1,5669 | 1,1889 |
| 1500 | 11,3821 | 11,1716 | 5,1348 | 4,6310 |
| 2000 | 26,6769 | 24,0502 | 11,9912 | 14,4759 |
| 2500 | 51,7488 | 53,1444 | 23,2098 | 23,5398 |
| 3000 | 89,0138 | 88,2979 | 39,8643 | 36,3688 |
Рис. 7.10. Зависимость ускорения от размера матриц при выполнении параллельного алгоритма Кэннона
Рис. 7.11. График зависимости экспериментального и теоретического времени выполнения алгоритма Кэннона на четырех процессорах
Умножение матриц при ленточной схеме разделения данных
Рассмотрим два параллельных алгоритма умножения матриц, в которых матрицы A и B разбиваются на непрерывные последовательности строк или столбцов (полосы).
Выделение информационных зависимостей
Для вычисления одной строки матрицы С необходимо, чтобы в каждой подзадаче содержалась строка матрицы А и был обеспечен доступ ко всем столбцам матрицы B. Возможные способы организации параллельных вычислений состоят в следующем.
1. Первый алгоритм. Алгоритм представляет собой итерационную процедуру, количество итераций которой совпадает с числом подзадач. На каждой итерации алгоритма каждая подзадача содержит по одной строке матрицы А и одному столбцу матрицы В. При выполнении итерации проводится скалярное умножение содержащихся в подзадачах строк и столбцов, что приводит к получению соответствующих элементов результирующей матрицы С. По завершении вычислений в конце каждой итерации столбцы матрицы В должны быть переданы между подзадачами с тем, чтобы в каждой подзадаче оказались новые столбцы матрицы В и могли быть вычислены новые элементы матрицы C. При этом данная передача столбцов между подзадачами должна быть организована таким образом, чтобы после завершения итераций алгоритма в каждой подзадаче последовательно оказались все столбцы матрицы В.
Возможная простая схема организации необходимой последовательности передач столбцов матрицы В между подзадачами состоит в представлении топологии информационных связей подзадач в виде кольцевой структуры. В этом случае на каждой итерации подзадача i, 0i<n, будет передавать свой столбец матрицы В подзадаче с номером i+1 (в соответствии с кольцевой структурой подзадача n-1 передает свои данные подзадаче с номером 0) – см. рис. 7.2. После выполнения всех итераций алгоритма необходимое условие будет обеспечено – в каждой подзадаче поочередно окажутся все столбцы матрицы В.
На рис. 7.2 представлены итерации алгоритма матричного умножения для случая, когда матрицы состоят из четырех строк и четырех столбцов (n=4). В начале вычислений в каждой подзадаче i, 0i<n, располагаются i-я строка матрицы A и i-й столбец матрицы B. В результате их перемножения подзадача получает элемент cii
результирующей матрицы С. Далее подзадачи осуществляют обмен столбцами, в ходе которого каждая подзадача передает свой столбец матрицы B следующей подзадаче в соответствии с кольцевой структурой информационных взаимодействий.
Итак, за основу параллельных вычислений для матричного умножения при блочном разделении данных принят подход, при котором базовые подзадачи отвечают за вычисления отдельных блоков матрицы C и при этом в подзадачах на каждой итерации расчетов располагается только по одному блоку исходных матриц A и B. Для нумерации подзадач будем использовать индексы размещаемых в подзадачах блоков матрицы C, т.е. подзадача (i,j) отвечает за вычисление блока Cij – тем самым, набор подзадач образует квадратную решетку, соответствующую структуре блочного представления матрицы C.
Возможный способ организации вычислений при таких условиях состоит в применении широко известного алгоритма Фокса (Fox) — см., например, [[34], [51]].
В соответствии с алгоритмом Фокса в ходе вычислений на каждой базовой подзадаче (i,j) располагается четыре матричных блока:
блок Cij матрицы C, вычисляемый подзадачей;блок Aij матрицы A, размещаемый в подзадаче перед началом вычислений;блоки A'ij , B'ij матриц A и B, получаемые подзадачей в ходе выполнения вычислений.
Выполнение параллельного метода включает:
этап инициализации, на котором каждой подзадаче (i,j) передаются блоки Aij, Bij и обнуляются блоки Cij на всех подзадачах;этап вычислений, в рамках которого на каждой итерации l, 0l<q, осуществляются следующие операции:для каждой строки i, 0i<q, блок Aij подзадачи (i,j) пересылается на все подзадачи той же строки i решетки; индекс j, определяющий положение подзадачи в строке, вычисляется в соответствии с выражением
где mod есть операция получения остатка от целочисленного деления; полученные в результаты пересылок блоки A'ij, B'ij каждой подзадачи (i,j) перемножаются и прибавляются к блоку Cijблоки B'ij каждой подзадачи (i,j) пересылаются подзадачам, являющимся соседями сверху в столбцах решетки подзадач (блоки подзадач из первой строки решетки пересылаются подзадачам последней строки решетки).
Рис. 7.6. Состояние блоков в каждой подзадаче в ходе выполнения итераций алгоритма Фокса
Для пояснения этих правил параллельного метода на рис. 7.6
приведено состояние блоков в каждой подзадаче в ходе выполнения итераций этапа вычислений (для решетки подзадач 2?2).
Отличие алгоритма Кэннона от метода Фокса состоит в изменении схемы начального распределения блоков перемножаемых матриц между подзадачами вычислительной системы. Начальное расположение блоков в алгоритме Кэннона подбирается таким образом, чтобы блоки в подзадачах могли бы быть перемножены без каких-либо дополнительных передач данных. При этом подобное распределение блоков может быть организовано таким образом, что перемещение блоков между подзадачами в ходе вычислений может осуществляться с использованием более простых коммуникационных операций.
С учетом высказанных замечаний этап инициализации алгоритма Кэннона включает выполнение следующих операций передач данных:
в каждую подзадачу (i,j) передаются блоки Aij, Bij;для каждой строки i решетки подзадач блоки матрицы A сдвигаются на (i-1) позиций влево;для каждого столбца j решетки подзадач блоки матрицы B сдвигаются на (j-1) позиций вверх.
Выполняемые при перераспределении матричных блоков процедуры передачи данных являются примером операции циклического сдвига – см. лекцию 3. Для пояснения используемого способа начального распределения данных на рис. 7.9 показан пример расположения блоков для решетки подзадач 3і3.
увеличить изображение
Рис. 7.9. Перераспределение блоков исходных матриц между процессорами при выполнении алгоритма Кэннона
В результате такого начального распределения в каждой базовой подзадаче будут располагаться блоки, которые могут быть перемножены без дополнительных операций передачи данных. Кроме того, получение всех последующих блоков для подзадач может быть обеспечено при помощи простых коммуникационных действий — после выполнения операции блочного умножения каждый блок матрицы A должен быть передан предшествующей подзадаче влево по строкам решетки подзадач, а каждый блок матрицы В – предшествующей подзадаче вверх по столбцам решетки. Как можно показать, последовательность таких циклических сдвигов и умножение получаемых блоков исходных матриц A и B приведет к получению в базовых подзадачах соответствующих блоков результирующей матрицы C.
Выполните реализацию двух ленточных алгоритмов
Выполните реализацию двух ленточных алгоритмов умножения матриц. Сравните время выполнения этих алгоритмов. Выполните реализацию алгоритма Кэннона. Постройте теоретические оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты реальных экспериментов с ранее полученными теоретическими оценками. Выполните реализацию блочных алгоритмов умножения матриц, которые могли бы быть выполнены для прямоугольных процессорных решеток общего вида.Выполните реализацию матричного умножения с использованием ранее разработанных программ умножения матрицы на вектор.
Алгоритм Гаусса
Метод Гаусса – широко известный прямой алгоритм решения систем линейных уравнений, для которых матрицы коэффициентов являются плотными. Если система линейных уравнений невырожденна, то метод Гаусса гарантирует нахождение решения с погрешностью, определяемой точностью машинных вычислений. Основная идея метода состоит в приведении матрицы А посредством эквивалентных преобразований (не меняющих решение системы (8.2)) к треугольному виду, после чего значения искомых неизвестных могут быть получены непосредственно в явном виде.
В подразделе дается общая характеристика метода Гаусса, достаточная для начального понимания алгоритма и позволяющая рассмотреть возможные способы параллельных вычислений при решении систем линейных уравнений. Более полное изложение алгоритма со строгим обсуждением вопросов точности получаемых решений может быть получено, например, в работах [[6], [22], [47]] и др.
Анализ эффективности
Оценим трудоемкость рассмотренного параллельного варианта метода Гаусса. Пусть, как и ранее, n есть порядок решаемой системы линейных уравнений, а p, p<n, обозначает число используемых процессоров. Тем самым, матрица коэффициентов А имеет размер n?n, и, соответственно, n/p есть размер полосы матрицы А на каждом процессоре.
Прежде всего, несложно показать, что общее время выполнения последовательного варианта метода Гаусса составляет:
|
(8.3) |
Определим теперь сложность параллельного варианта метода Гаусса. При выполнении прямого хода алгоритма на каждой итерации для выбора ведущей строки каждый процессор должен осуществить выбор максимального значения в столбце с исключаемой неизвестной в пределах своей полосы. Начальный размер полос на процессорах равен n/p, но по мере исключения неизвестных количество строк в полосах для обработки постепенно сокращается. Текущий размер полос приближенно можно оценить как (n-i)/p, где i, 0i<n-1, есть номер выполняемой итерации прямого хода метода Гаусса. Далее после сбора полученных максимальных значений, определения и рассылки ведущей строки каждый процессор должен выполнить вычитание ведущей строки из каждой строки оставшейся части строк своей полосы матрицы A. Количество элементов строки, подлежащих обработке, также сокращается при исключении неизвестных, и текущее число элементов строки для вычислений оценивается величиной (n-i). Тем самым, сложность процедуры вычитания строк оценивается как 2·(n-i) операций (перед вычитанием ведущая строка умножается на масштабирующую величину aik/aii). С учетом выполняемого количества итераций общее число операций параллельного варианта прямого хода метода Гаусса определяется выражением:
На каждой итерации обратного хода алгоритма Гаусса после рассылки вычисленного значения очередной неизвестной каждый процессор должен обновить значения правых частей для всех строк, расположенных на этом процессоре. Отсюда следует, что трудоемкость параллельного варианта обратного хода алгоритма Гаусса оценивается как величина:
Выберем для дальнейшего анализа эффективности получаемых параллельных вычислений параллельный алгоритм матрично-векторного умножения при ленточном горизонтальном разделении матрицы и при полном дублировании всех обрабатываемых векторов.
Трудоемкость последовательного метода сопряженных градиентов была уже определена ранее в (8.12).
Определим время выполнения параллельной реализации метода сопряженных градиентов. Вычислительная сложность параллельных операций умножения матрицы на вектор при использовании схемы ленточного горизонтального разделения матрицы составляет:
(см. лекцию 6).
Как результат, при условии дублирования всех вычислений над векторами общая вычислительная сложность параллельного варианта метода сопряженных градиентов равна:
С учетом полученных оценок показатели ускорения и эффективности параллельного алгоритма могут быть выражены при помощи соотношений:
Рассмотрев построенные показатели, можно отметить, что балансировка вычислительной нагрузки между процессорами в целом является достаточно равномерной.
Уточним теперь приведенные выражения – учтем длительность выполняемых вычислительных операций и оценим трудоемкость операции передачи данных между процессорами. Как можно заметить, информационное взаимодействие процессоров возникает только при выполнении операций умножения матрицы на вектор. С учетом результатов лекции 6 коммуникационная сложность рассматриваемых параллельных вычислений равна:
где – латентность, ? – пропускная способность сети передачи данных, а w есть размер элемента упорядочиваемых данных в байтах.
Окончательно, время выполнения параллельного варианта метода сопряженных градиентов для решения систем линейных уравнений принимает вид:
|
(8.13) |
где ? есть время выполнения базовой вычислительной операции.
Контрольные вопросы
Что представляет собой система линейных уравнений? Какие типы систем вам известны? Какие методы могут быть использованы для решения систем разных типов?В чем состоит постановка задачи решения системы линейных уравнений?В чем идея параллельной реализации метода Гаусса?Какие информационные взаимодействия имеются между базовыми подзадачами для параллельного варианта метода Гаусса?Каковы показатели эффективности для параллельного варианта метода Гаусса?В чем состоит схема программной реализации параллельного варианта метода Гаусса?В чем состоит идея параллельной реализации метода сопряженных градиентов?Какой из алгоритмов обладает большей коммуникационной сложностью?
Краткий обзор лекции
Данная лекция посвящена проблеме параллельных вычислений при решении систем линейных уравнений. Изложение учебного материала проводится с использованием двух широко известных алгоритмов: метода Гаусса как примера прямого алгоритма решения задачи и итерационного метода сопряженных градиентов.
Параллельный вариант метода Гаусса (подраздел 8.2) основывается на ленточном разделении матрицы между процессорами с использованием циклической схемы распределения строк, что позволяет сбалансировать вычислительную нагрузку. Для определения параллельного варианта метода проводится полный цикл проектирования – определяются базовые подзадачи, выделяются информационные взаимодействия, обсуждаются вопросы масштабирования, выводятся оценки показателей эффективности, предлагается схема программной реализации и приводятся результаты вычислительных экспериментов. Как видно из графика, представленного на рис. 8.8, параллельный алгоритм Гаусса демонстрирует достаточно высокие показатели ускорения и эффективности.
Важный момент при рассмотрении параллельного варианта метода сопряженных градиентов (подраздел 8.3) состоит в том, что способ параллельных вычислений для этого метода может быть получен через параллельные алгоритмы выполняемых вычислительных действий – операций умножения матрицы на вектор, скалярного произведения векторов, сложения и вычитания векторов. Выбранный для учебного примера подход состоит в распараллеливании вычислительно наиболее трудоемкой операции умножения матрицы на вектор, в то время как все вычисления над векторами дублируются на каждом процессоре для уменьшения числа выполняемых операций передачи данных — показатели эффективности подобного способа организации параллельных вычислений представлены на рис. 8.8.
Рис. 8.8. Ускорение параллельных алгоритмов решения системы линейных уравнений с размером матрицы 3000x3000
Масштабирование и распределение подзадач по процессорам
Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью и сбалансированным объемом передаваемых данных. В случае когда размер матрицы, описывающей систему линейных уравнений, оказывается большим, чем число доступных процессоров (т.е. p<n), базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько строк матрицы. Однако применение использованной в лекциях 6 и 7
последовательной схемы разделения данных для параллельного решения систем линейных уравнений приведет к неравномерной вычислительной нагрузке процессоров: по мере исключения (на прямом ходе) или определения (на обратном ходе) неизвестных в методе Гаусса для все большей части процессоров все необходимые вычисления будут завершены и процессоры окажутся простаивающими. Возможное решение проблемы балансировки вычислений может состоять в использовании ленточной циклической схемы (см. лекцию 6) для распределения данных между укрупненными подзадачами. В этом случае матрица A делится на наборы (полосы) строк вида (см. рис. 8.2):
Рис. 8.2. Пример использования ленточной циклической схемы разделения строк матрицы между тремя процессорами
Сопоставив схему разделения данных и порядок выполнения вычислений в методе Гаусса, можно отметить, что использование циклического способа формирования полос позволяет обеспечить лучшую балансировку вычислительной нагрузки между подзадачами.
Распределение подзадач между процессорами должно учитывать характер выполняемых в методе Гаусса коммуникационных операций. Основным видом информационного взаимодействия подзадач является операция передачи данных от одного процессора всем процессорам вычислительной системы. Как результат, для эффективной реализации требуемых информационных взаимодействий между базовыми подзадачами топология сети передачи данных должны иметь структуру гиперкуба или полного графа.
Метод сопряженных градиентов
Рассмотрим теперь совершенно иной подход к решению систем линейных уравнений, при котором к искомому точному решению x*
системы Ax=b строится последовательность приближенных решений x0, x1, ..., xk, ... При этом процесс вычислений организуется таким способом, что каждое очередное приближение дает оценку точного решения со все уменьшающейся погрешностью, и при продолжении расчетов оценка точного решения может быть получена с любой требуемой точностью. Подобные итерационные методы решения систем линейных уравнений широко используются в практике вычислений. К преимуществам итерационных методов можно отнести меньший объем (по сравнению, например, с методом Гаусса) необходимых вычислений для решения разреженных систем линейных уравнений, возможность более быстрого получения начальных приближений искомого решения, наличие эффективных способов организации параллельных вычислений. Дополнительная информация с описанием таких методов, рассмотрением вопросов сходимости и точности получаемых решений может быть получена, например, в [[6], [22]].
Рис. 8.4. График зависимости экспериментального и теоретического времени проведения эксперимента на двух процессорах от объема исходных данных
Одним из наиболее известных итерационных алгоритмов является метод сопряженных градиентов, который может быть применен для решения симметричной положительно определенной системы линейных уравнений большой размерности.
Напомним, что матрица А является симметричной, если она совпадает со своей транспонированной матрицей, т.е. А=АТ. Матрица А называется положительно определенной, если для любого вектора x справедливо: xTAx>0.
Известно (см., например, [[6], [22]]), что после выполнения n итераций алгоритма (n есть порядок решаемой системы линейных уравнений), очередное приближение xn совпадает с точным решением.
Обратный ход алгоритма Гаусса
После приведения матрицы коэффициентов к верхнему треугольному виду становится возможным определение значений неизвестных. Из последнего уравнения преобразованной системы может быть вычислено значение переменной xn-1, после этого из предпоследнего уравнения становится возможным определение переменной xn-2 и т.д. В общем виде выполняемые вычисления при обратном ходе метода Гаусса могут быть представлены при помощи соотношений:
Поясним, как и ранее, выполнение обратного хода метода Гаусса на примере рассмотренной в п. 8.2.1.1 системы линейных уравнений
Из последнего уравнения системы можно определить, что неизвестная x2 имеет значение 3. В результате становится возможным разрешение второго уравнения и определение значение неизвестной x1=13, т.е.
На последней итерации обратного хода метода Гаусса определяется значение неизвестной x0, равное -44.
С учетом последующего параллельного выполнения можно отметить, что вычисление получаемых значений неизвестных может выполняться сразу во всех уравнениях системы (и эти действия могут выполняться в уравнениях одновременно и независимо друг от друга). Так, в рассматриваемом примере после определения значения неизвестной x2 система уравнений может быть приведена к виду
Вычислительная сложность обратного хода алгоритма Гаусса составляет O(n2).
Обзор литературы
Проблема численного решения систем линейных уравнений широко рассматривается в литературе. Для учебного рассмотрения могут быть рекомендованы работы [[6], [13], [51], [63]]. Широкое обсуждение вопросов параллельных вычислений для решения данного типа задач выполнено в работах [[22], [30]].
При рассмотрении вопросов программной реализации параллельных методов может быть рекомендована работа [[23]]. В ней рассматривается хорошо известная и широко используемая в практике параллельных вычислений программная библиотека численных методов ScaLAPACK.
Определение подзадач
При внимательном рассмотрении метода Гаусса можно заметить, что все вычисления сводятся к однотипным вычислительным операциям над строками матрицы коэффициентов системы линейных уравнений. Как результат, в основу параллельной реализации алгоритма Гаусса может быть положен принцип распараллеливания по данным. В качестве базовой подзадачи можно принять тогда все вычисления, связанные с обработкой одной строки матрицы A и соответствующего элемента вектора b.
Организация параллельных вычислений
При разработке параллельного варианта метода сопряженных градиентов для решения систем линейных уравнений в первую очередь следует учесть, что выполнение итераций метода осуществляется последовательно и, тем самым, наиболее целесообразный подход состоит в распараллеливании вычислений, реализуемых в ходе выполнения итераций.
Анализ соотношений (8.8) – (8.11) показывает, что основные вычисления, выполняемые в соответствии с методом, состоят в умножении матрицы A на векторы x и d, и, как результат, при организации параллельных вычислений может быть полностью использован материал, изложенный в лекции 6.
Дополнительные вычисления в (8.8) – (8.11), имеющие меньший порядок сложности, представляют собой различные операции обработки векторов (скалярное произведение, сложение и вычитание, умножение на скаляр). Организация таких вычислений, конечно же, должна быть согласована с выбранным параллельным способом выполнения операции умножения матрицы на вектор. Общие же рекомендации могут состоять в следующем: при малом размере векторов можно применить дублирование векторов между процессорами, при большом порядке решаемой системы линейных уравнений более целесообразно осуществлять блочное разделение векторов.
Последовательный алгоритм
Если матрица A симметричная и положительно определенная, то функция
|
(8.6) |
имеет единственный минимум, который достигается в точке x*, совпадающей с решением системы линейных уравнений (8.2). Метод сопряженных градиентов является одним из широкого класса итерационных алгоритмов, которые позволяют найти решение (8.2) путем минимизации функции q(x).
Итерация метода сопряженных градиентов состоит в вычислении очередного приближения к точному решению в соответствии с правилом:
|
(8.7) |
Тем самым, новое значение приближения xk вычисляется с учетом приближения, построенного на предыдущем шаге xk-1, скалярного шага sk и вектора направления dk.
Перед выполнением первой итерации векторы x0 и d0 полагаются равными нулю, а для вектора g0 устанавливается значение, равное - b. Далее каждая итерация для вычисления очередного значения приближения xk включает выполнение четырех шагов:
Шаг 1: Вычисление градиента:
|
(8.8) |
Шаг 2: Вычисление вектора направления:
где (gT,g) есть скалярное произведение векторов.
Шаг 3: Вычисление величины смещения по выбранному направлению:
|
(8.10) |
Шаг 4: Вычисление нового приближения:
|
(8.11) |
Как можно заметить, данные выражения включают две операции умножения матрицы на вектор, четыре операции скалярного произведения и пять операций над векторами. Как результат, общее количество операций, выполняемых на одной итерации, составляет
t1=2n2+13n.
Как уже отмечалось ранее, для нахождения точного решения системы линейных уравнений с положительно определенной симметричной матрицей необходимо выполнить n итераций. Таким образом, для нахождения решения системы необходимо выполнить
|
(8.12) |
и, тем самым, сложность алгоритма имеет порядок O(n3).
Поясним выполнение метода сопряженных градиентов на примере решения системы линейных уравнений вида:
3x0 -x1 = 3, -x0 +3x1 = 7.
Эта система уравнений второго порядка обладает симметричной положительно определенной матрицей, для нахождения точного решения этой системы достаточно провести всего две итерации метода.
На первой итерации было получено значение градиента g1=(-3,-7), значение вектора направления d1=(3, 7), значение величины смещения s1=0,439. Соответственно, очередное приближение к точному решению системы x1=(1,318, 3,076).
На второй итерации было получено значение градиента g2=(-2,121, 0,909), значение вектора направления d2=(2,397, -0,266), а величина смещения – s2=0,284. Очередное приближение совпадает с точным решением системы x2=(2, 3).
На рис. 8.5 представлена последовательность приближений к точному решению, построенная методом сопряженных градиентов (в качестве начального приближения x0 выбрана точка (0,0)).
Рис. 8.5. Итерации метода сопряженных градиентов при решении системы второго порядка
Метод Гаусса основывается на возможности выполнения преобразований линейных уравнений, которые не меняют при этом решения рассматриваемой системы (такие преобразования носят наименование эквивалентных). К числу таких преобразований относятся:
умножение любого из уравнений на ненулевую константу;перестановка уравнений;прибавление к уравнению любого другого уравнения системы.
Метод Гаусса включает последовательное выполнение двух этапов. На первом этапе – прямой ход метода Гаусса – исходная система линейных уравнений при помощи последовательного исключения неизвестных приводится к верхнему треугольному виду
Ux=c,
где матрица коэффициентов получаемой системы имеет вид
На обратном ходе метода Гаусса (второй этап алгоритма) осуществляется определение значений неизвестных. Из последнего уравнения преобразованной системы может быть вычислено значение переменной xn-1, после этого из предпоследнего уравнения становится возможным определение переменной xn-2 и т.д.
Постановка задачи
Линейное уравнение с n неизвестными x0, x1, ѕ, xn-1 может быть определено при помощи выражения
|
(8.1) |
где величины a0,a1,...,an-1 и b представляют собой постоянные значения.
Множество из n линейных уравнений
|
(8.2) |
называется системой линейных уравнений или линейной системой. В более кратком (матричном) виде система может представлена как
Ax = b,
где A=(ai,j) есть вещественная матрица размера n?n, а векторы b и x состоят из n элементов.
Под задачей решения системы линейных уравнений для заданных матрицы А и вектора b обычно понимается нахождение значения вектора неизвестных x, при котором выполняются все уравнения системы.
Алгоритм Гаусса решения систем линейных
|
// Программа 8.1. — Алгоритм Гаусса решения систем линейных уравнений int ProcNum; // Число доступных процессоров int ProcRank; // Ранг текущего процессора int *pParallelPivotPos; // Номера строк, которые были выбраны ведущими int *pProcPivotIter; // Номера итераций, на которых строки // процессора использовались в качестве ведущих void main(int argc, char* argv[]) { double* pMatrix; // Матрица линейной системы double* pVector; // Вектор правых частей линейной системы double* pResult; // Вектор неизвестных double *pProcRows; // Строки матрицы A double *pProcVector; // Блок вектора b double *pProcResult; // Блок вектора x int Size; // Размер матрицы и векторов int RowNum; // Количество строк матрицы double start, finish, duration; setvbuf(stdout, 0, _IONBF, 0); MPI_Init ( &argc, &argv ); MPI_Comm_rank ( MPI_COMM_WORLD, &ProcRank ); MPI_Comm_size ( MPI_COMM_WORLD, &ProcNum ); if (ProcRank == 0) printf("Параллельный метод Гаусса для решения систем линейных уравнений\n"); // Выделение памяти и инициализация данных ProcessInitialization(pMatrix, pVector, pResult, pProcRows, pProcVector, pProcResult, Size, RowNum); // Распределение исходных данных DataDistribution(pMatrix, pProcRows, pVector, pProcVector, Size, RowNum); // Выполнение параллельного алгоритма Гаусса ParallelResultCalculation(pProcRows, pProcVector, pProcResult, Size, RowNum); // Сбор найденного вектора неизвестных на ведущем процессе ResultCollection(pProcResult, pResult); // Завершение процесса вычислений ProcessTermination(pMatrix, pVector, pResult, pProcRows, pProcVector, pProcResult); MPI_Finalize(); } |
| Пример 8.1. |
| Закрыть окно |
// Программа 8.1. — Алгоритм Гаусса решения систем линейных уравнений
int ProcNum; // Число доступных процессоров
int ProcRank; // Ранг текущего процессора
int *pParallelPivotPos; // Номера строк, которые были выбраны ведущими
int *pProcPivotIter; // Номера итераций, на которых строки
// процессора использовались в качестве ведущих
void main(int argc, char* argv[]) {
double* pMatrix; // Матрица линейной системы
double* pVector; // Вектор правых частей линейной системы
double* pResult; // Вектор неизвестных
double *pProcRows; // Строки матрицы A
double *pProcVector; // Блок вектора b
double *pProcResult; // Блок вектора x
int Size; // Размер матрицы и векторов
int RowNum; // Количество строк матрицы
double start, finish, duration;
setvbuf(stdout, 0, _IONBF, 0);
MPI_Init ( &argc, &argv );
MPI_Comm_rank ( MPI_COMM_WORLD, &ProcRank );
MPI_Comm_size ( MPI_COMM_WORLD, &ProcNum );
if (ProcRank == 0)
printf("Параллельный метод Гаусса для решения систем линейных уравнений\n");
// Выделение памяти и инициализация данных
ProcessInitialization(pMatrix, pVector, pResult,
pProcRows, pProcVector, pProcResult, Size, RowNum);
// Распределение исходных данных
DataDistribution(pMatrix, pProcRows, pVector, pProcVector,
Size, RowNum);
// Выполнение параллельного алгоритма Гаусса
ParallelResultCalculation(pProcRows, pProcVector, pProcResult, Size,
RowNum);
// Сбор найденного вектора неизвестных на ведущем процессе
ResultCollection(pProcResult, pResult);
// Завершение процесса вычислений
ProcessTermination(pMatrix, pVector, pResult, pProcRows,
pProcVector, pProcResult);
MPI_Finalize();
}
Прямой ход алгоритма Гаусса
Прямой ход метода Гаусса состоит в последовательном исключении неизвестных в уравнениях решаемой системы линейных уравнений. На итерации i, 0i<n-1, метода производится исключение неизвестной i для всех уравнений с номерами k, большими i (т.е. i<kn-1). Для этого из этих уравнений осуществляется вычитание строки i, умноженной на константу (aki/aii), с тем чтобы результирующий коэффициент при неизвестной xi в строках оказался нулевым – все необходимые вычисления могут быть определены при помощи соотношений:
(следует отметить, что аналогичные вычисления выполняются и над вектором b).
Поясним выполнение прямого хода метода Гаусса на примере системы линейных уравнений вида:
На первой итерации производится исключение неизвестной x0 из второй и третьей строки. Для этого из этих строк нужно вычесть первую строку, умноженную соответственно на 2 и 1. После этих преобразований система уравнений принимает вид:
В результате остается выполнить последнюю итерацию и исключить неизвестную x1 из третьего уравнения. Для этого необходимо вычесть вторую строку, и в окончательной форме система имеет следующий вид:
На рис. 8.1 представлена общая схема состояния данных на i-й итерации прямого хода алгоритма Гаусса. Все коэффициенты при неизвестных, расположенные ниже главной диагонали и левее столбца i, уже являются нулевыми. На i-й итерации прямого хода метода Гаусса осуществляется обнуление коэффициентов столбца i, расположенных ниже главной диагонали, путем вычитания строки i, умноженной на нужную ненулевую константу. После проведения (n-1) подобной итерации матрица, определяющая систему линейных уравнений, становится приведенной к верхнему треугольному виду.
Рис. 8.1. Итерация прямого хода алгоритма Гаусса
При выполнении прямого хода метода Гаусса строка, которая используется для исключения неизвестных, носит наименование ведущей, а диагональный элемент ведущей строки называется ведущим элементом. Как можно заметить, выполнение вычислений является возможным только, если ведущий элемент имеет ненулевое значение. Более того, если ведущий элемент ai,i имеет малое значение, то деление и умножение строк на этот элемент может приводить к накоплению вычислительной погрешности и вычислительной неустойчивости алгоритма.
Возможный способ избежать подобной проблемы может состоять в следующем: при выполнении каждой очередной итерации прямого хода метода Гаусса следует определить коэффициент с максимальным значением по абсолютной величине в столбце, соответствующем исключаемой неизвестной, т.е.
и выбрать в качестве ведущей строку, в которой этот коэффициент располагается (данная схема выбора ведущего значения носит наименование метода главных элементов).
Вычислительная сложность прямого хода алгоритма Гаусса с выбором ведущей строки имеет порядок O(n3).
Программная реализация
Рассмотрим возможный вариант параллельной реализации метода Гаусса для решения систем линейных уравнений. Следует отметить, что для получения более простого вида программы для разделения матрицы между процессами используется ленточная последовательная схема (полосы матрицы образуют последовательные наборы соседних строк матрицы). В описании программы реализация отдельных модулей не приводится, если их отсутствие не оказывает влияния на понимание общей схемы параллельных вычислений.
1. Главная функция программы. Реализует логику работы алгоритма, последовательно вызывает необходимые подпрограммы.
Пример 8.1.
(html, txt)
Следует пояснить использование дополнительных массивов. Элементы массива pParallelPivotPos определяют номера строк матрицы, выбираемых в качестве ведущих, по итерациям прямого хода метода Гаусса. Именно в этом порядке должны выполняться затем итерации обратного хода для определения значений неизвестных системы линейных уравнений. Массив pParallelPivotPos является глобальным, и любое его изменение в одном из процессов требует выполнения операции рассылки измененных данных всем остальным процессам программы.
Элементы массива pProcPivotIter определяют номера итераций прямого хода метода Гаусса, на которых строки процесса использовались в качестве ведущих (т.е. строка i процесса выбиралась ведущей на итерации pProcPivotIter[i]). Начальное значение элементов массива устанавливается нулевым, и, тем самым, нулевое значение элемента массива pProcPivotIter[i] является признаком того, что строка i процесса все еще подлежит обработке. Кроме того, важно отметить, что запоминаемые в элементах массива pProcPivotIter номера итераций дополнительно означают и номера неизвестных, для определения которых будут использованы соответствующие строки уравнения. Массив pProcPivotIter является локальным для каждого процесса.
Функция ProcessInitialization определяет исходные данные решаемой задачи (число неизвестных), выделяет память для хранения данных, осуществляет ввод матрицы коэффициентов системы линейных уравнений и вектора правых частей (или формирует эти данные при помощи какого-либо датчика случайных чисел).
Функция DataDistribution реализует распределение матрицы линейной системы и вектора правых частей между процессорами вычислительной системы.
Функция ResultCollection осуществляет сбор со всех процессов отдельных частей вектора неизвестных.
Функция ProcessTermination выполняет необходимый вывод результатов решения задачи и освобождает всю ранее выделенную память для хранения данных.
Реализация всех перечисленных функций может быть выполнена по аналогии с ранее рассмотренными примерами и предоставляется читателю в качестве самостоятельного упражнения.
2. Функция ParallelResultCalculation. Реализует логику работы параллельного алгоритма Гаусса, последовательно вызывает функции, выполняющие прямой и обратный ход метода.
// Функция для параллельного выполнения метода Гаусса void ParallelResultCalculation (double* pProcRows, double* pProcVector, double* pProcResult, int Size, int RowNum) { ParallelGaussianElimination (pProcRows, pProcVector, Size, RowNum); ParallelBackSubstitution (pProcRows, pProcVector, pProcResult, Size, RowNum); }
3. Функция ParallelGaussianElimination. Функция выполняет параллельный вариант прямого хода алгоритма Гаусса.
Пример 8.2.
(html, txt)
Функция ParallelEliminateColumns проводит вычитание ведущей строки из строк процесса, которые еще не использовались в качестве ведущих (т.е. для которых элементы массива pProcPivotIter равны нулю).
4. Функция ParallelBackSubstitution. Функция реализует параллельный вариант обратного хода Гаусса.
Пример 8.3.
(html, txt)
Функция FindBackPivotRow определяет строку, из которой можно вычислить значение очередного элемента результирующего вектора. Номер этой строки хранится в массиве pParallelPivotIter. По номеру функция FindBackPivotRow определяет номер процесса, на котором эта строка хранится, и номер этой строки в полосе pProcRows этого процесса.
Решение систем линейных уравнений
Системы линейных уравнений возникают при решении ряда прикладных задач, описываемых дифференциальными, интегральными или системами нелинейных (трансцендентных) уравнений. Они могут появляться также в задачах математического программирования, статистической обработки данных, аппроксимации функций, при дискретизации краевых дифференциальных задач методом конечных разностей или методом конечных элементов и др.
Матрицы коэффициентов систем линейных уравнений могут иметь различные структуру и свойства. Матрицы решаемых систем могут быть плотными, и их порядок может достигать несколько тысяч строк и столбцов. При решении многих задач могут появляться системы, обладающие симметричными положительно определенными ленточными матрицами с порядком в десятки тысяч и шириной ленты в несколько тысяч элементов. И, наконец, при рассмотрении большого ряда задач могут возникать системы линейных уравнений с разреженными матрицами с порядком в миллионы строк и столбцов.
Результаты вычислительных экспериментов
Вычислительные эксперименты для оценки эффективности параллельного варианта метода Гаусса для решения систем линейных уравнений проводились при условиях, указанных в п. 6.5.5, и состоят в следующем.
Эксперименты производились на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server (см. п. 1.2.3).
Для оценки длительности ? базовой скалярной операции проводилось решение системы линейных уравнений при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины ? было получено значение 4,7 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности a и пропускной способности b соответственно 47 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.
Результаты вычислительных экспериментов приведены в таблице 8.1. Эксперименты выполнялись с использованием двух, четырех и восьми процессоров.
| 500 | 0,36 | 0,3302 | 1,0901 | 0,5170 | 0,6963 | 0,7504 | 0,4796 |
| 1000 | 3,313 | 1,5950 | 2,0770 | 1,6152 | 2,0511 | 1,8715 | 1,7701 |
| 1500 | 11,437 | 4,1788 | 2,7368 | 3,8802 | 2,9474 | 3,7567 | 3,0443 |
| 2000 | 26,688 | 9,3432 | 2,8563 | 7,2590 | 3,6765 | 7,3713 | 3,6204 |
| 2500 | 50,125 | 16,9860 | 2,9509 | 11,9957 | 4,1785 | 11,6530 | 4,3014 |
| 3000 | 85,485 | 28,4948 | 3,0000 | 19,1255 | 4,4696 | 17,6864 | 4,8333 |
Рис. 8.3. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма Гаусса для разных размеров систем линейных уравнений
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (8.5) приведено в таблице 8.2 и на рис. 8.4.
Вычислительные эксперименты для оценки эффективности параллельного варианта метода сопряженных градиентов для решения систем линейных уравнений проводились при условиях, указанных в п. 8.2.7.
Результаты вычислительных экспериментов приведены в таблице 8.3. Эксперименты проводились на вычислительных системах, состоящих из двух, четырех и восьми процессоров.
| 500 | 0,5 | 0,4634 | 1,0787 | 0,4706 | 1,0623 | 1,3020 | 0,3840 |
| 1000 | 8,14 | 3,9207 | 2,0761 | 3,6354 | 2,2390 | 3,5092 | 2,3195 |
| 1500 | 31,391 | 17,9505 | 1,7487 | 14,4102 | 2,1783 | 20,2001 | 1,5539 |
| 2000 | 92,36 | 51,3204 | 1,7996 | 40,7451 | 2,2667 | 37,9319 | 2,4348 |
| 2500 | 170,549 | 125,3005 | 1,3611 | 85,0761 | 2,0046 | 87,2626 | 1,9544 |
| 3000 | 363,476 | 223,3364 | 1,6274 | 146,1308 | 2,4873 | 134,1359 | 2,7097 |
Рис. 8.6. Зависимость ускорения от количества процессоров при выполнении параллельного метода сопряженных градиентов для решения систем линейных уравнений
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (8.13) приведено в таблице 8.4 и на рис. 8.7.
| 500 | 1,3042 | 0,4634 | 0,6607 | 0,4706 | 0,3390 | 1,3020 |
| 1000 | 10,3713 | 3,9207 | 5,2194 | 3,6354 | 2,6435 | 3,5092 |
| 34,9333 | 17,9505 | 17,5424 | 14,4102 | 8,8470 | 20,2001 | |
| 2000 | 82,7220 | 51,3204 | 41,4954 | 40,7451 | 20,8822 | 37,9319 |
| 2500 | 161,4695 | 125,3005 | 80,9446 | 85,0761 | 40,6823 | 87,2626 |
| 3000 | 278,9077 | 223,3364 | 139,7560 | 146,1308 | 70,1801 | 134,1359 |
Рис. 8.7. График зависимости экспериментального и теоретического времени проведения эксперимента на четырех процессорах от объема исходных данных
Выделение информационных зависимостей
Рассмотрим общую схему параллельных вычислений и возникающие при этом информационные зависимости между базовыми подзадачами.
Для выполнения прямого хода метода Гаусса необходимо осуществить (n-1) итерацию по исключению неизвестных для преобразования матрицы коэффициентов A к верхнему треугольному виду.
Выполнение итерации i, 0i<n-1, прямого хода метода Гаусса включает ряд последовательных действий. Прежде всего, в самом начале итерации необходимо выбрать ведущую строку, которая при использовании метода главных элементов определяется поиском строки с наибольшим по абсолютной величине значением среди элементов столбца i, соответствующего исключаемой переменной xi. Поскольку строки матрицы A распределены по подзадачам, для поиска максимального значения подзадачи с номерами k, k<i, должны обменяться своими элементами при исключаемой переменной xi. После сбора всех необходимых данных в каждой подзадаче может быть определено, какая из подзадач содержит ведущую строку и какое значение является ведущим элементом.
Далее для продолжения вычислений ведущая подзадача должна разослать свою строку матрицы A и соответствующий элемент вектора b всем остальным подзадачам с номерами k, k<i. Получив ведущую строку, подзадачи выполняют вычитание строк, обеспечивая тем самым исключение соответствующей неизвестной xi.
При выполнении обратного хода метода Гаусса подзадачи выполняют необходимые вычисления для нахождения значения неизвестных. Как только какая-либо подзадача i, 0i<n-1, определяет значение своей переменной xi, это значение должно быть разослано всем подзадачам с номерами k, k<i. Далее подзадачи подставляют полученное значение новой неизвестной и выполняют корректировку значений для элементов вектора b.
Выполните анализ эффективности параллельных вычислений
Выполните анализ эффективности параллельных вычислений в отдельности для прямого и обратного этапов метода Гаусса. Оцените, на каком этапе происходит большее снижение показателей.Выполните разработку параллельного варианта метода Гаусса при вертикальном разбиении матрицы по столбцам. Постройте теоретические оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты выполненных экспериментов с ранее полученными теоретическими оценками.Выполните реализацию параллельного метода сопряженных градиентов. Постройте теоретические оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты выполненных экспериментов с ранее полученными теоретическими оценками.Выполните разработку параллельных вариантов методов Якоби и Зейделя решения систем линейных уравнений (см., например, Kumar (1994)). Постройте теоретические оценки времени работы этого алгоритма с учетом параметров используемой вычислительной системы. Проведите вычислительные эксперименты. Сравните результаты выполненных экспериментов с ранее полученными теоретическими оценками.
