Алгоритм чет-нечетной перестановки
Алгоритм пузырьковой сортировки в прямом виде достаточно сложен для распараллеливания – сравнение пар значений упорядочиваемого набора данных происходит строго последовательно. В связи с этим для параллельного применения используется не сам этот алгоритм, а его модификация, известная в литературе как метод чет-нечетной перестановки (the odd-even transposition method) – см., например, [[51]]. Суть модификации состоит в том, что в алгоритм сортировки вводятся два разных правила выполнения итераций метода: в зависимости от четности или нечетности номера итерации сортировки для обработки выбираются элементы с четными или нечетными индексами соответственно, сравнение выделяемых значений всегда осуществляется с их правыми соседними элементами. Таким образом, на всех нечетных итерациях сравниваются пары
(a1, a2), (a3, a4), ..., (an-1,an) (при четном n),
а на четных итерациях обрабатываются элементы
(a2, a3), (a4, a5), ..., (an-2,an-1).
После n-кратного повторения итераций сортировки исходный набор данных оказывается упорядоченным.
Алгоритм 9.2. Последовательный алгоритм чет-нечетной перестановки
// Алгоритм 9.2 // Последовательный алгоритм чет-нечетной перестановки void OddEvenSort(double A[], int n) { for (int i = 1; i < n; i++) { if (i % 2 == 1) { // нечетная итерация for (int j = 0; j < n/2 - 2; j++) compare_exchange(A[2*j + 1], A[2*j + 2]); if (n % 2 == 1) // сравнение последней пары при нечетном n compare_exchange(A[n - 2], A[n - 1]); } else // четная итерация for (int j = 1; j < n/2 - 1; j++) compare_exchange(A[2*j], A[2*j + 1]); } }
Анализ эффективности
При анализе эффективности, как и ранее, вначале проведем общую оценку сложности рассмотренного параллельного алгоритма сортировки, а затем дополним полученные соотношения показателями трудоемкости выполняемых коммуникационных операций.
Определим первоначально трудоемкость последовательных вычислений. При рассмотрении данного вопроса алгоритм пузырьковой сортировки позволяет продемонстрировать следующий важный момент. Как уже отмечалось в начале этой лекции, использованный для распараллеливания последовательный метод упорядочивания данных характеризуется квадратичной зависимостью сложности от числа упорядочиваемых данных, т.е. T1~n2. Однако применение подобной оценки сложности последовательного алгоритма приведет к искажению исходного целевого назначения критериев качества параллельных вычислений – показатели эффективности в этом случае будут характеризовать используемый способ параллельного выполнения данного конкретного метода сортировки, а не результативность использования параллелизма для задачи упорядочивания данных в целом как таковой. Различие состоит в том, что для сортировки могут быть применены более эффективные последовательные алгоритмы, трудоемкость которых имеет порядок
|
(9.1) |
и чтобы сравнить, насколько быстрее могут быть упорядочены данные при использовании параллельных вычислений, в обязательном порядке должна применяться именно эта оценка сложности. Как основной результат приведенных рассуждений, можно сформулировать утверждение о том, что при определении показателей ускорения и эффективности параллельных вычислений в качестве оценки сложности последовательного способа решения рассматриваемой задачи следует использовать трудоемкость наилучших последовательных алгоритмов. Параллельные методы решения задач должны сравниваться с наиболее быстродействующими последовательными способами вычислений!
Определим теперь сложность рассмотренного параллельного алгоритма упорядочивания данных. Как отмечалось ранее, на начальной стадии работы метода каждый процессор проводит упорядочивание своих блоков данных (размер блоков при равномерном распределении данных равен n/p).
Для оценки эффективности параллельного аналога алгоритма Шелла могут быть использованы соотношения, полученные для параллельного метода пузырьковой сортировки (см. п. 9.3.5). При этом следует только учесть двухэтапность алгоритма Шелла – с учетом данной особенности общее время выполнения нового параллельного метода может быть определено при помощи выражения
|
(9.7) |
Как можно заметить, эффективность параллельного варианта сортировки Шелла существенно зависит от значения L – при малом значении величины L новый параллельный способ сортировки выполняется быстрее, чем ранее рассмотренный алгоритм чет-нечетной перестановки.
Оценим трудоемкость рассмотренного параллельного метода. Пусть у нас имеется N-мерный гиперкуб, состоящий из p=2N процессоров, где p<n.
Эффективность параллельного метода быстрой сортировки, как и в последовательном варианте, во многом зависит от правильности выбора значений ведущих элементов. Определение общего правила для выбора этих значений представляется затруднительным. Сложность такого выбора может быть снижена, если выполнить упорядочение локальных блоков процессоров перед началом сортировки и обеспечить однородное распределение сортируемых данных между процессорами вычислительной системы.
Определим вначале вычислительную сложность алгоритма сортировки. На каждой из log2p итераций сортировки каждый процессор осуществляет деление блока относительно ведущего элемента, сложность этой операции составляет n/p операций (будем предполагать, что на каждой итерации сортировки каждый блок делится на равные по размеру части).
При завершении вычислений процессор выполняет сортировку своих блоков, что может быть выполнено при использовании быстрых алгоритмов за (n/p)log2(n/p) операций.
Таким образом, общее время вычислений параллельного алгоритма быстрой сортировки составляет
|
(9.8) |
где ? есть время выполнения базовой операции перестановки.
Рассмотрим теперь сложность выполняемых коммуникационных операций. Общее количество межпроцессорных обменов для рассылки ведущего элемента на N-мерном гиперкубе может быть ограничено оценкой
|
(9.9) |
При используемых предположениях (выбор ведущих элементов осуществляется наилучшим образом) количество итераций алгоритма равно log2p, а объем передаваемых данных между процессорами всегда равен половине блока, т.е. величине (n/p)/2. При таких условиях коммуникационная сложность параллельного алгоритма быстрой сортировки определяется при помощи соотношения:
|
(9.10) |
где – латентность, ? – пропускная способность сети, а w есть размер элемента набора в байтах.
С учетом всех полученных соотношений общая трудоемкость алгоритма оказывается равной
|
(9.11) |
Оценим трудоемкость рассмотренного параллельного метода. Пусть, как и ранее, n есть количество сортируемых данных, p, p<n, обозначает число используемых процессоров и, соответственно, n/p есть размер блоков данных на процессорах.
В течение первого этапа алгоритма каждый процессор сортирует свой блок данных с помощью быстрой сортировки, тем самым, длительность выполняемых при этом операций равна
|
(9.13) |
где ? есть время выполнения базовой операции сортировки.
На втором этапе алгоритма один из процессоров собирает наборы из p элементов со всех остальных процессоров, выполняет слияние всех полученных данных (общее количество элементов составляет p2), формирует набор из p-1 ведущих элементов и рассылает полученный набор всем остальным процессорам. С учетом всех перечисленных действий общая длительность второго этапа составляет
|
(9.14) |
(в приведенном соотношении выделенные подвыражения соответствуют четырем перечисленным действиям алгоритма); здесь, как и ранее, – латентность, ? – пропускная способность сети передачи данных, а w есть размер элемента упорядочиваемых данных в байтах.
В ходе выполнения третьего этапа алгоритма каждый процессор разделяет свои элементы относительно ведущих элементов на p частей (общее количество операций для этого может быть ограничено величиной n/p). Далее все процессоры выполняют рассылку сформированных частей блоков между собой – оценка трудоемкости такой коммуникационной операции рассмотрена в лекции 3 при представлении топологии вычислительной сети в виде гиперкуба. Как было показано, выполнение такой операции может быть осуществлено за log2p шагов, на каждом из которых каждый процессор передает и получает сообщение из (n/p)/2 элементов. Как результат, общая трудоемкость третьего этапа алгоритма может быть оценена как
|
(9.15) |
На четвертом этапе алгоритма каждый процессор выполняет слияние p отсортированных частей в один объединенный блок. Оценка трудоемкости такой операции уже проводилась при рассмотрении второго этапа, и, тем самым, длительность выполнения процедуры слияния составляет
|
(9.16) |
С учетом всех полученных соотношений общее время выполнения алгоритма сортировки с использованием регулярного набора образцов составляет
|
(9.17) |
Контрольные вопросы
В чем состоит постановка задачи сортировки данных?Приведите несколько примеров алгоритмов сортировки. Какова вычислительная сложность приведенных алгоритмов?Какая операция является базовой для задачи сортировки данных?В чем суть параллельного обобщения базовой операции задачи сортировки данных?Что представляет собой алгоритм чет-нечетной перестановки?В чем состоит параллельный вариант алгоритма Шелла? Каковы основные отличия этого параллельного алгоритма сортировки от метода чет-нечетной перестановки?Что представляет собой параллельный вариант алгоритма быстрой сортировки?Что зависит от правильного выбора ведущего элемента для параллельного алгоритма быстрой сортировки? Какие способы выбора ведущего элемента могут быть предложены?Для каких топологий могут применяться рассмотренные алгоритмы сортировки?В чем состоит алгоритм сортировки с использованием регулярного набора образцов?
Краткий обзор лекции
В лекции рассматривается часто встречающаяся в приложениях задача упорядочения данных, для решения которой в рамках данного учебного материала выбраны широко известные алгоритмы пузырьковой сортировки, сортировки Шелла и быстрой сортировки. При изложении методов сортировки основное внимание уделяется возможным способам распараллеливания алгоритмов, анализу эффективности и сравнению получаемых теоретических оценок с результатами выполненных вычислительных экспериментов.
Алгоритм пузырьковой сортировки (подраздел 9.3) в исходном виде практически не поддается распараллеливанию в силу последовательного выполнения основных итераций метода. Для введения необходимого параллелизма рассматривается обобщенный вариант алгоритма – метод чет-нечетной перестановки. Суть обобщения состоит в том, что в алгоритм сортировки вводятся два разных правила выполнения итераций метода в зависимости от четности номера итерации сортировки. Сравнения пар значений упорядочиваемого набора данных на итерациях метода чет-нечетной перестановки являются независимыми и могут быть выполнены параллельно.
Для алгоритма Шелла (подраздел 9.4) рассматривается схема распараллеливания при представлении топологии сети в виде гиперкуба. При таком представлении топологии оказывается возможной организация взаимодействия процессоров, расположенных далеко друг от друга при линейной нумерации. Как правило, такая организация вычислений позволяет уменьшить количество выполняемых итераций алгоритма сортировки.
Для алгоритма быстрой сортировки (подраздел 9.5) приводятся три схемы распараллеливания. Первые две схемы также основываются на представлении топологии сети в виде гиперкуба. Основная итерация вычислений состоит в выборе одним из процессоров ведущего элемента, который далее рассылается всем остальным процессорам гиперкуба. После получения ведущего элемента процессоры проводят разделение своих блоков, и получаемые части блоков передаются между попарно связанными процессорами. В результате выполнения подобной итерации исходный гиперкуб оказывается разделенным на 2 гиперкуба меньшей размерности, к которым, в свою очередь, может быть применена приведенная выше схема вычислений.
При применении алгоритма быстрой сортировки одним из основных моментов является правильность выбора ведущего элемента. Оптимальная стратегия состоит в выборе такого значения ведущего элемента, при котором данные на процессорах разделяются на части одинакового размера. В общем случае при произвольно сгенерированных исходных данных достижение такой ситуации является достаточно сложной задачей. В первой схеме предлагается выбирать ведущий элемент, например, как среднее арифметическое элементов на управляющем процессоре. Во второй схеме (обобщенный алгоритм быстрой сортировки) данные на каждом процессоре предварительно упорядочиваются с тем, чтобы взять средний элемент блока данных как ведущее значение. Третья схема (сортировка с использованием регулярного набора образцов) распараллеливания алгоритма быстрой сортировки основывается на многоуровневой схеме формирования множества ведущих элементов. Такой подход может быть применен для произвольного количества процессоров, позволяет избежать множества пересылок данных и приводит, как правило, к более равномерному распределению данных между процессорами.
Масштабирование и распределение подзадач по процессорам
Как отмечалось ранее, количество подзадач соответствует числу имеющихся процессоров, и поэтому необходимости в проведении масштабирования вычислений не возникает. Исходное распределение блоков упорядочиваемого набора данных по процессорам может быть выбрано совершенно произвольным образом. Для эффективного выполнения рассмотренного параллельного алгоритма сортировки нужно, чтобы процессоры с соседними номерами имели прямые линии связи.
Масштабирование параллельных вычислений
Рассмотренное параллельное обобщение базовой операции сортировки может быть надлежащим образом адаптировано и для случая p<n, когда количество процессоров меньше числа упорядочиваемых значений. В данной ситуации каждый процессор будет содержать уже не единственное значение, а часть (блок размера n/p) сортируемого набора данных.
Определим в качестве результата выполнения параллельного алгоритма сортировки такое состояние упорядочиваемого набора данных, при котором имеющиеся на процессорах данные упорядочены, а порядок распределения блоков по процессорам соответствует линейному порядку нумерации (т. е. значение последнего элемента на процессоре Pi меньше значения первого элемента на процессоре Pi+1, где 0i<p-1 или равно ему).
Блоки обычно упорядочиваются в самом начале сортировки на каждом процессоре в отдельности при помощи какого-либо быстрого алгоритма (начальная стадия параллельной сортировки). Далее, следуя схеме одноэлементного сравнения, взаимодействие пары процессоров Pi и Pi+1 для совместного упорядочения содержимого блоков Ai и Ai+1 может быть осуществлено следующим образом:
выполнить взаимообмен блоков между процессорами Pi и Pi+1;объединить блоки Ai и Ai+1 на каждом процессоре в один отсортированный блок двойного размера (при исходной упорядоченности блоков Ai и Ai+1 процедура их объединения сводится к быстрой операции слияния упорядоченных наборов данных);разделить полученный двойной блок на две равные части и оставить одну из этих частей (например, с меньшими значениями данных) на процессоре Pi, а другую часть (с большими значениями, соответственно) – на процессоре Pi+1
Рассмотренная процедура обычно именуется в литературе операцией "сравнить и разделить" (compare-split). Следует отметить, что сформированные в результате такой процедуры блоки на процессорах Pi и Pi+1 совпадают по размеру с исходными блоками Ai и Ai+1 и все значения, расположенные на процессоре Pi, не превышают значений на процессоре Pi+1.
Определенная выше операция "сравнить и разделить" может быть использована в качестве базовой подзадачи для организации параллельных вычислений. Как следует из построения, количество таких подзадач параметрически зависит от числа имеющихся процессоров, и, таким образом, проблема масштабирования вычислений для параллельных алгоритмов сортировки практически отсутствует. Вместе с тем следует отметить, что относящиеся к подзадачам блоки данных изменяются в ходе выполнения сортировки. В простых случаях размер блоков данных в подзадачах остается неизмененным. В более сложных ситуациях (как, например, в алгоритме быстрой сортировки – см. подраздел 9.5) объем располагаемых на процессорах данных может различаться, что может приводить к нарушению равномерной вычислительной загрузки процессоров.
Обобщенный алгоритм быстрой сортировки
В обобщенном алгоритме быстрой сортировки (the HyperQuickSort algorithm) в дополнение к обычному методу быстрой сортировки предлагается конкретный способ выбора ведущих элементов. Суть предложения состоит в том, что сортировка располагаемых на процессорах блоков происходит в самом начале выполнения вычислений. Кроме того, для поддержки упорядоченности в ходе вычислений процессоры должны выполнять операции слияния частей блоков, получаемых после разделения. Как результат, в силу упорядоченности блоков, при выполнении алгоритма быстрой сортировки в качестве ведущего элемента целесообразнее будет выбирать средний элемент какого-либо блока (например, на первом процессоре вычислительной системы). Выбираемый подобным образом ведущий элемент в отдельных случаях может оказаться более близким к реальному среднему значению всего сортируемого набора, чем какое-либо другое произвольно выбранное значение.
Все остальные действия в новом рассматриваемом алгоритме выполняются в соответствии с обычным методом быстрой сортировки. Более подробное описание данного способа распараллеливания быстрой сортировки может быть получено, например, в [[63]].
При анализе эффективности обобщенного алгоритма можно воспользоваться соотношением (9.11). Следует только учесть, что на каждой итерации метода теперь выполняется операция слияния частей блоков (будем, как и ранее, предполагать, что их размер одинаков и равен (n/p)/2). Кроме того, в силу упорядоченности блоков может быть усовершенствована процедура деления – вместо перебора всех элементов блока теперь достаточно будет выполнить для ведущего элемента бинарный поиск в блоке. С учетом всех высказанных замечаний трудоемкость обобщенного алгоритма быстрой сортировки может быть выражена при помощи следующего выражения:
|
(9.12) |
Обзор литературы
Возможные способы решения задачи упорядочения данных широко обсуждаются в литературе; один из наиболее полных обзоров алгоритмов сортировки содержится в работе [[50]], среди последних изданий может быть рекомендована работа [[26]].
Параллельные варианты алгоритма пузырьковой сортировки и сортировки Шелла рассматриваются в [[51]].
Схемы распараллеливания быстрой сортировки при представлении топологии сети передачи данных в виде гиперкуба описаны в [[51], [63]]. Сортировка с использованием регулярного набора образцов представлена в работе [[63]].
Полезной при рассмотрении вопросов параллельных вычислений для сортировки данных может оказаться работа [[17]].
Определение подзадач и выделение информационных зависимостей
Получение параллельного варианта для метода чет-нечетной перестановки уже не представляет каких-либо затруднений – сравнения пар значений на итерациях сортировки являются независимыми и могут быть выполнены параллельно. В случае p<n, когда количество процессоров меньше числа упорядочиваемых значений, процессоры содержат блоки данных размера n/p и в качестве базовой подзадачи может быть использована операция "сравнить и разделить" (см. подраздел 9.2).
Алгоритм 9.3. Параллельный алгоритм чет-нечетной перестановки
// Алгоритм 9.3 // Параллельный алгоритм чет-нечетной перестановки ParallelOddEvenSort(double A[], int n) { int id = GetProcId(); // номер процесса int np = GetProcNum(); // количество процессов for (int i = 0; i < np; i++ ) { if (i % 2 == 1) { // нечетная итерация if (id % 2 == 1) { // нечетный номер процесса // Cравнение-обмен с процессом — соседом справа if (id < np - 1) compare_split_min(id + 1); } else // Cравнение-обмен с процессом — соседом слева if (id > 0) compare_split_max(id - 1); } else { // четная итерация if(id % 2 == 0) { // четный номер процесса // Cравнение-обмен с процессом — соседом справа if (id < np - 1) compare_split_min(id + 1); } else // Cравнение-обмен с процессом — соседом слева compare_split_max(id - 1); } } }
Для пояснения такого параллельного способа сортировки в табл. 9.1 приведен пример упорядочения данных при n=16, p=4 (т.е. блок значений на каждом процессоре содержит n/p=4 элемента). В первом столбце таблицы приводится номер и тип итерации метода, перечисляются пары процессов, для которых параллельно выполняются операции "сравнить и разделить". Взаимодействующие пары процессов выделены в таблице рамкой. Для каждого шага сортировки показано состояние упорядочиваемого набора данных до и после выполнения итерации.
В общем случае выполнение параллельного метода может быть прекращено, если в течение каких-либо двух последовательных итераций сортировки состояние упорядочиваемого набора данных не было изменено.
Как результат, общее количество итераций может быть сокращено, и для фиксации таких моментов необходимо введение некоторого управляющего процессора, который определял бы состояние набора данных после выполнения каждой итерации сортировки. Однако трудоемкость такой коммуникационной операции (сборка на одном процессоре сообщений от всех процессоров) может оказаться столь значительной, что весь эффект от возможного сокращения итераций сортировки будет поглощен затратами на реализацию операций межпроцессорной передачи данных.
| Исходные данные | 13 55 59 88 | 29 43 71 85 | 2 18 40 75 | 4 14 22 43 |
| 1 нечет (1, 2), (3, 4) | 13 55 59 88 | 29 43 71 85 | 2 18 40 75 | 4 14 22 43 |
| 13 29 43 55 | 59 71 85 88 | 2 4 14 18 | 22 40 43 75 | |
| 2 чет (2, 3) | 13 29 43 55 | 59 71 85 88 | 2 4 14 18 | 22 40 43 75 |
| 13 29 43 55 | 2 4 14 18 | 59 71 85 88 | 22 40 43 75 | |
| 3 нечет (1, 2), (3, 4) | 13 29 43 55 | 2 4 14 18 | 59 71 85 88 | 22 40 43 75 |
| 2 4 13 14 | 18 29 43 55 | 22 40 43 59 | 71 75 85 88 | |
| 4 чет (2, 3) | 2 4 13 14 | 18 29 43 55 | 22 40 43 59 | 71 75 85 88 |
| 2 4 13 14 | 18 22 29 40 | 43 43 55 59 | 71 75 85 88 |
Организация параллельных вычислений
Для алгоритма Шелла может быть предложен параллельный аналог метода (см., например, [[51]]), если топология коммуникационной сети может быть эффективно представлена в виде N-мерного гиперкуба (т. е. количество процессоров равно p = 2N). Выполнение сортировки в таком случае может быть разделено на два последовательных этапа. На первом этапе (N итераций) осуществляется взаимодействие процессоров, являющихся соседними в структуре гиперкуба (но эти процессоры могут оказаться далекими при линейной нумерации; для установления соответствия двух систем нумерации процессоров может быть использован, как и ранее, код Грея – см. лекцию 3). Формирование пар процессоров, взаимодействующих между собой при выполнении операции "сравнить и разделить", может быть обеспечено при помощи следующего простого правила – на каждой итерации i, 0i<N, парными становятся процессоры, у которых различие в битовых представлениях их номеров имеется только в позиции N-i.
Второй этап состоит в реализации обычных итераций параллельного алгоритма чет-нечетной перестановки. Итерации данного этапа выполняются до прекращения фактического изменения сортируемого набора, и, тем самым, общее количество L таких итераций может быть различным – от 2 до p.
На рис. 9.3 показан пример сортировки массива из 16 элементов с помощью рассмотренного способа (процессоры показаны кружками, номера процессоров даны в битовом представлении). Нужно заметить, что данные оказываются упорядоченными уже после первого этапа и нет необходимости выполнять чет-нечетную перестановку.
Рис. 9.3. Пример работы алгоритма Шелла для 4 процессоров
С учетом представленного описания параллельного варианта алгоритма Шелла базовая подзадача для организации параллельных вычислений, как и ранее (см. п. 9.3.3), может быть определена на основе операции "сравнить и разделить". Как результат, количество подзадач всегда совпадает с числом имеющихся процессоров (размер блоков данных в подзадачах равен n/p) и не возникает проблемы масштабирования. Распределение блоков упорядочиваемого набора данных по процессорам должно быть выбрано с учетом возможности эффективного выполнения операций "сравнить и разделить" при представлении топологии сети передачи данных в виде гиперкуба.
Алгоритм сортировки с использованием регулярного набора образцов (the parallel sorting by regular sampling) также является обобщением метода быстрой сортировки (см., например, в [[63]]).
Упорядочивание данных в соответствии с данным вариантом алгоритма быстрой сортировки осуществляется в ходе выполнения следующих четырех этапов:
на первом этапе сортировки производится упорядочивание имеющихся на процессорах блоков. Данная операция может быть выполнена каждым процессором независимо друг от друга при помощи обычного алгоритма быстрой сортировки; далее каждый процессор формирует набор из элементов своих блоков с индексами 0, m, 2m, ...,(p-1)m, где m=n/p2;на втором этапе выполнения алгоритма все сформированные на процессорах наборы данных собираются на одном из процессоров системы и объединяются в ходе последовательного слияния в одно упорядоченное множество. Далее из полученного множества значений из элементов с индексами
формируется новый набор ведущих элементов, который передается всем используемым процессорам. В завершение этапа каждый процессор выполняет разделение своего блока на p частей с использованием полученного набора ведущих значений;на третьем этапе сортировки каждый процессор осуществляет рассылку выделенных ранее частей своего блока всем остальным процессорам системы; рассылка выполняется в соответствии с порядком нумерации – часть j, 0j<p, каждого блока пересылается процессору с номером j;на четвертом этапе выполнения алгоритма каждый процессор выполняет слияние p полученных частей в один отсортированный блок.
По завершении четвертого этапа исходный набор данных становится отсортированным.
На рис. 9.11 приведен пример сортировки массива данных с помощью алгоритма, описанного выше. Следует отметить, что число процессоров для данного алгоритма может быть произвольным, в данном примере оно равно 3.
Рис. 9.11. Пример работы алгоритма сортировки с использованием регулярного набора образцов
Параллельное обобщение алгоритма быстрой сортировки (см., например, [63]) наиболее простым способом может быть получено, если топология коммуникационной сети может быть эффективно представлена в виде N-мерного гиперкуба (т.е. p=2N). Пусть, как и ранее, исходный набор данных распределен между процессорами блоками одинакового размера n/p; результирующее расположение блоков должно соответствовать нумерации процессоров гиперкуба. Возможный способ выполнения первой итерации параллельного метода при таких условиях может состоять в следующем:
выбрать каким-либо образом ведущий элемент и разослать его по всем процессорам системы (например, в качестве ведущего элемента можно взять среднее арифметическое элементов, расположенных на выбранном ведущем процессоре);разделить на каждом процессоре имеющийся блок данных на две части с использованием полученного ведущего элемента;образовать пары процессоров, для которых битовое представление номеров отличается только в позиции N, и осуществить взаимообмен данными между этими процессорами.
В результате выполнения такой итерации сортировки исходный набор оказывается разделенным на две части, одна из которых (со значениями меньшими, чем значение ведущего элемента) располагается на процессорах, в битовом представлении номеров которых бит N равен 0. Таких процессоров всего p/2, и, таким образом, исходный N-мерный гиперкуб также оказывается разделенным на два гиперкуба размерности N-1. К этим подгиперкубам, в свою очередь, может быть параллельно применена описанная выше процедура. После N-кратного повторения подобных итераций для завершения сортировки достаточно упорядочить блоки данных, получившиеся на каждом отдельном процессоре вычислительной системы.
Для пояснения на рис. 9.6 представлен пример упорядочивания данных при n=16, p=4 (т.е. блок каждого процессора содержит 4 элемента). На этом рисунке процессоры изображены в виде прямоугольников, внутри которых показано содержимое упорядочиваемых блоков данных; значения блоков приводятся в начале и при завершении каждой итерации сортировки. Взаимодействующие пары процессоров соединены двунаправленными стрелками. Для разделения данных выбирались наилучшие значения ведущих элементов: на первой итерации для всех процессоров использовалось значение 0, на второй итерации для пары процессоров 0, 1 ведущий элемент равен -5, для пары процессоров 2, 3 это значение было принято равным 4.
Рис. 9.6. Пример упорядочивания данных параллельным методом быстрой сортировки (без результатов локальной сортировки блоков)
Как и ранее, в качестве базовой подзадачи для организации параллельных вычислений может быть выбрана операция "сравнить и разделить", а количество подзадач совпадает с числом используемых процессоров. Распределение подзадач по процессорам должно производиться с учетом возможности эффективного выполнения алгоритма при представлении топологии сети передачи данных в виде гиперкуба.
Параллельные методы сортировки
Сортировка является одной из типовых проблем обработки данных и обычно понимается как задача размещения элементов неупорядоченного набора значений
в порядке монотонного возрастания или убывания
(здесь и далее все пояснения для краткости будут даваться только на примере упорядочивания данных по возрастанию).
Возможные способы решения этой задачи широко обсуждаются в литературе; один из наиболее полных обзоров алгоритмов сортировки содержится в работе [[50]], среди последних изданий может быть рекомендована работа [[26]].
Вычислительная трудоемкость процедуры упорядочивания является достаточно высокой. Так, для ряда известных простых методов (пузырьковая сортировка, сортировка включением и др.) количество необходимых операций определяется квадратичной зависимостью от числа упорядочиваемых данных
Для более эффективных алгоритмов (сортировка слиянием, сортировка Шелла, быстрая сортировка) трудоемкость определяется величиной T \sim n \log_2 n.
Данное выражение дает также нижнюю оценку необходимого количества операций для упорядочивания набора из n значений; алгоритмы с меньшей трудоемкостью могут быть получены только для частных вариантов задачи.
Ускорение сортировки может быть обеспечено при использовании нескольких (p>1) процессоров. Исходный упорядочиваемый набор в этом случае разделяется между процессорами; в ходе сортировки данные пересылаются между процессорами и сравниваются между собой. Результирующий (упорядоченный) набор, как правило, также разделен между процессорами; при этом для систематизации такого разделения для процессоров вводится та или иная система последовательной нумерации и обычно требуется, чтобы при завершении сортировки значения, располагаемые на процессорах с меньшими номерами, не превышали значений процессоров с большими номерами.
Оставляя подробный анализ проблемы сортировки для отдельного рассмотрения, здесь основное внимание мы уделим изучению параллельных способов выполнения для ряда широко известных методов внутренней сортировки, когда все упорядочиваемые данные могут быть размещены полностью в оперативной памяти ЭВМ.
Последовательный алгоритм
Последовательный алгоритм пузырьковой сортировки (the bubble sort algorithm) (см., например, [[26], [50]]) сравнивает и обменивает соседние элементы в последовательности, которую нужно отсортировать. Для последовательности
(a1, a2, ..., an)
алгоритм сначала выполняет n-1 базовых операций "сравнения-обмена" для последовательных пар элементов
(a1, a2), (a2, a3), ..., (an-1, an).
В результате после первой итерации алгоритма самый большой элемент перемещается ("всплывает") в конец последовательности. Далее последний элемент в преобразованной последовательности может быть исключен из рассмотрения, и описанная выше процедура применяется к оставшейся части последовательности
(a'1, a'2, ..., a'n-1).
Как можно увидеть, последовательность будет отсортирована после n-1 итерации. Эффективность пузырьковой сортировки может быть улучшена, если завершать алгоритм в случае отсутствия каких- либо изменений сортируемой последовательности данных в ходе какой-либо итерации сортировки.
Алгоритм 9.1. Последовательный алгоритм пузырьковой сортировки
// Алгоритм 9.1. // Последовательный алгоритм пузырьковой сортировки void BubbleSort(double A[], int n) { for (int i = 0; i < n - 1; i++) for (int j = 0; j < n - i; j++) compare_exchange(A[j], A[j + 1]); }
Общая идея сортировки Шелла (the Shell sort) (см., например, [26, 50]) состоит в сравнении на начальных стадиях сортировки пар значений, располагаемых достаточно далеко друг от друга в упорядочиваемом наборе данных. Такая модификация метода сортировки позволяет быстро переставлять далекие неупорядоченные пары значений (сортировка таких пар обычно требует большого количества перестановок, если используется сравнение только соседних элементов).
Общая схема метода состоит в следующем. На первом шаге алгоритма происходит упорядочивание элементов n/2 пар (ai, an/2+i) для 1in/2. Далее, на втором шаге упорядочиваются элементы в n/4 группах из четырех элементов (ai, an/4+i, an/2+i, a3n/4+i) для 1in/4. На третьем шаге упорядочиваются элементы уже в n/4 группах из восьми элементов и т.д. На последнем шаге упорядочиваются элементы сразу во всем массиве (a1, a2, ..., an). На каждом шаге для упорядочивания элементов в группах применяется метод сортировки вставками. Как можно заметить, общее количество итераций алгоритма Шелла равно log2n.
В более полном виде алгоритм Шелла может быть представлен следующим образом.
Алгоритм 9.4. Последовательный алгоритм сортировки Шелла
// Алгоритм 9.4 // Последовательный алгоритм сортировки Шелла void ShellSort(double A[], int n) { int incr = n / 2; while (incr > 0) { for (int i = incr + 1; i < n; i++) { int j = i - incr; while (j > 0) if (A[j] > A[j + incr]) { swap(A[j], A[j + incr]); j = j - incr; } else j = 0; } incr = incr/2; } }
При общем рассмотрении алгоритма быстрой сортировки (the quick sort algorithm), предложенной Хоаром (C.A.R. Hoare), прежде всего следует отметить, что этот метод основывается на последовательном разделении сортируемого набора данных на блоки меньшего размера таким образом, что между значениями разных блоков обеспечивается отношение упорядоченности (для любой пары блоков все значения одного из этих блоков не превышают значений другого блока). На первой итерации метода осуществляется деление исходного набора данных на первые две части – для организации такого деления выбирается некоторый ведущий элемент и все значения набора, меньшие ведущего элемента, переносятся в первый формируемый блок, все остальные значения образуют второй блок набора. На второй итерации сортировки описанные правила применяются рекурсивно для обоих сформированных блоков и т.д. При надлежащем выборе ведущих элементов после выполнения log2n итераций исходный массив данных оказывается упорядоченным. Более подробное изложение метода может быть получено, например, в [[26], [50]].
Эффективность быстрой сортировки в значительной степени определяется правильностью выбора ведущих элементов при формировании блоков. В худшем случае трудоемкость метода имеет тот же порядок сложности, что и пузырьковая сортировка (т.е. T1~n2). При оптимальном выборе ведущих элементов, когда разделение каждого блока происходит на равные по размеру части, трудоемкость алгоритма совпадает с быстродействием наиболее эффективных способов сортировки . В среднем случае количество операций, выполняемых алгоритмом быстрой сортировки, определяется выражением (см., например, [[26], [50]]):
Общая схема алгоритма быстрой сортировки может быть представлена в следующем виде (в качестве ведущего элемента выбирается первый элемент упорядочиваемого набора данных).
Алгоритм 9.5. Последовательный алгоритм быстрой сортировки
// Алгоритм 9.5 // Последовательный алгоритм быстрой сортировки void QuickSort(double A[], int i1, int i2) { if (i1 < i2) { double pivot = A[i1]; int is = i1; for (int i = i1 + 1; i < i2; i++) if (A[i] Ј pivot) { is = is + 1; swap(A[is], A[i]); } swap(A[i1], A[is]); QuickSort(A, i1, is); QuickSort(A, is + 1, i2); } }
Функция для выполнения обощенного алгоритма
|
// Функция для выполнения обощенного алгоритма быстрой сортировки void ParallelHyperQuickSort ( double *pProcData, int ProcDataSize) { MPI_Status status; int CommProcRank; // Ранг процессора, с которым выполняется // взаимодействие double *pData, // Часть блока, остающаяся на процессоре *pSendData, // Часть блока, передаваемая процессору // CommProcRank *pRecvData, // Часть блока, получаемая от процессора // CommProcRank *pMergeData; // Блок данных, получаемый после слияния int DataSize, SendDataSize, RecvDataSize, MergeDataSize; int HypercubeDim = (int)ceil(log(ProcNum)/log(2)); // размерность гиперкуба int Mask = ProcNum; double Pivot; // Первоначальная сортировка блоков данных на каждом процессоре LocalDataSort(pProcData, ProcDataSize); // Итерации обобщенной быстрой сортировки for (int i = HypercubeDim; i > 0; i-- ) { // Определение ведущего значения и его рассылка всем процессорам PivotDistribution(pProcData, ProcDataSize, HypercubeDim, Mask, i,&Pivot); Mask = Mask >> 1; // Определение границы разделения блока int pos = GetProcDataDivisionPos(pProcData, ProcDataSize, Pivot); // Разделение блока на части if ( ( (rank & Mask) >> (i - 1) ) == 0 ) { // старший бит = 0 pSendData = &pProcData[pos + 1]; SendDataSize = ProcDataSize - pos – 1; if (SendDataSize < 0) SendDataSize = 0; CommProcRank = ProcRank + Mask pData = &pProcData[0]; DataSize = pos + 1; } else { // старший бит = 1 pSendData = &pProcData[0]; SendDataSize = pos + 1; if (SendDataSize > ProcDataSize) SendDataSize = pos; CommProcRank = ProcRank – Mask pData = &pProcData[pos + 1]; DataSize = ProcDataSize - pos - 1; if (DataSize < 0) DataSize = 0; } // Пересылка размеров частей блоков данных MPI_Sendrecv(&SendDataSize, 1, MPI_INT, CommProcRank, 0, &RecvDataSize, 1, MPI_INT, CommProcRank, 0, MPI_COMM_WORLD, &status); // Пересылка частей блоков данных pRecvData = new double[RecvDataSize]; MPI_Sendrecv(pSendData, SendDataSize, MPI_DOUBLE, CommProcRank, 0, pRecvData, RecvDataSize, MPI_DOUBLE, CommProcRank, 0, MPI_COMM_WORLD, &status); // Слияние частей MergeDataSize = DataSize + RecvDataSize; pMergeData = new double[MergeDataSize]; DataMerge(pMergeData, pMergeData, pData, DataSize, pRecvData, RecvDataSize); delete [] pProcData; delete [] pRecvData; pProcData = pMergeData; ProcDataSize = MergeDataSize; } } |
| Пример 9.1. |
| Закрыть окно |
// Функция для выполнения обощенного алгоритма быстрой сортировки
void ParallelHyperQuickSort ( double *pProcData,
int ProcDataSize) {
MPI_Status status;
int CommProcRank; // Ранг процессора, с которым выполняется
// взаимодействие
double *pData, // Часть блока, остающаяся на процессоре
*pSendData, // Часть блока, передаваемая процессору
// CommProcRank
*pRecvData, // Часть блока, получаемая от процессора
// CommProcRank
*pMergeData; // Блок данных, получаемый после слияния
int DataSize, SendDataSize, RecvDataSize, MergeDataSize;
int HypercubeDim = (int)ceil(log(ProcNum)/log(2));
// размерность гиперкуба
int Mask = ProcNum;
double Pivot;
// Первоначальная сортировка блоков данных на каждом процессоре
LocalDataSort(pProcData, ProcDataSize);
// Итерации обобщенной быстрой сортировки
for (int i = HypercubeDim; i > 0; i-- ) {
// Определение ведущего значения и его рассылка всем процессорам
PivotDistribution(pProcData, ProcDataSize, HypercubeDim,
Mask, i,&Pivot);
Mask = Mask >> 1;
// Определение границы разделения блока
int pos = GetProcDataDivisionPos(pProcData, ProcDataSize, Pivot);
// Разделение блока на части
if ( ( (rank & Mask) >> (i - 1) ) == 0 ) { // старший бит = 0
pSendData = &pProcData[pos + 1];
SendDataSize = ProcDataSize - pos – 1;
if (SendDataSize < 0) SendDataSize = 0;
CommProcRank = ProcRank + Mask
pData = &pProcData[0];
DataSize = pos + 1;
}
else { // старший бит = 1
pSendData = &pProcData[0];
SendDataSize = pos + 1;
if (SendDataSize > ProcDataSize) SendDataSize = pos;
CommProcRank = ProcRank – Mask
pData = &pProcData[pos + 1];
DataSize = ProcDataSize - pos - 1;
if (DataSize < 0) DataSize = 0;
}
// Пересылка размеров частей блоков данных
MPI_Sendrecv(&SendDataSize, 1, MPI_INT, CommProcRank, 0,
&RecvDataSize, 1, MPI_INT, CommProcRank, 0, MPI_COMM_WORLD, &status);
// Пересылка частей блоков данных
pRecvData = new double[RecvDataSize];
MPI_Sendrecv(pSendData, SendDataSize, MPI_DOUBLE,
CommProcRank, 0, pRecvData, RecvDataSize, MPI_DOUBLE,
CommProcRank, 0, MPI_COMM_WORLD, &status);
// Слияние частей
MergeDataSize = DataSize + RecvDataSize;
pMergeData = new double[MergeDataSize];
DataMerge(pMergeData, pMergeData, pData, DataSize,
pRecvData, RecvDataSize);
delete [] pProcData;
delete [] pRecvData;
pProcData = pMergeData;
ProcDataSize = MergeDataSize;
}
}
Принципы распараллеливания
При внимательном рассмотрении способов упорядочивания данных, применяемых в алгоритмах сортировки, можно заметить, что многие методы основаны на применении одной и той же базовой операции "сравнить и переставить" (compare-exchange), состоящей в сравнении той или иной пары значений из сортируемого набора данных и перестановке этих значений, если их порядок не соответствует условиям сортировки.
Пример 9.1. Операция "сравнить и переставить"
// Базовая операция "сравнить и переставить" if ( A[i] > A[j] ) { temp = A[i]; A[i] = A[j]; A[j] = temp; }
Целенаправленное применение данной операции позволяет упорядочить данные; в способах выбора пар значений для сравнения, собственно, и проявляется различие алгоритмов сортировки.
Для параллельного обобщения выделенной базовой операции сортировки рассмотрим первоначально ситуацию, когда количество процессоров совпадает с числом сортируемых значений (т. е. p=n) и на каждом из процессоров содержится только по одному значению исходного набора данных. Тогда сравнение значений ai и aj, располагаемых соответственно на процессорах Pi и Pj, можно организовать следующим образом (параллельное обобщение базовой операции сортировки):
выполнить взаимообмен имеющихся на процессорах Pi и Pj
значений (с сохранением на этих процессорах исходных элементов);сравнить на каждом процессоре Pi и Pj получившиеся одинаковые пары значений (ai, aj); результаты сравнения используются для разделения данных между процессорами – на одном процессоре (например, Pi) остается меньший элемент, другой процессор (т. е. Pj) запоминает для дальнейшей обработки большее значение пары
Представим возможный вариант параллельной программы обобщенной быстрой сортировки. При этом реализация отдельных модулей не приводится, если их отсутствие не оказывает влияния на понимание общей схемы параллельных вычислений.
1. Главная функция программы. Реализует логику работы алгоритма, последовательно вызывает необходимые подпрограммы.
// Программа 9.1. // Обобщенная быстрая сортировка int ProcRank; // Ранг текущего процесса int ProcNum; // Количество процессов int main(int argc, char *argv[]) { double *pProcData; // Блок данных процесса int ProcDataSize; // Размер блока данных
MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank); MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);
// Инициализация данных и их распределение между процессами ProcessInitialization(&pProcData, &ProcDataSize);
// Параллельная сортировка ParallelHyperQuickSort(pProcData, ProcDataSize);
// Завершение вычислений процесса ProcessTermination(pProcData, ProcDataSize);
MPI_Finalize(); }
Функция ProcessInitialization определяет исходные данные решаемой задачи (размер сортируемого массива), выделяет память для хранения данных, осуществляет генерацию сортируемого массива (например, при помощи датчика случайных чисел) и распределяет его между процессами.
Функция ProcessTermination осуществляет необходимый вывод результатов решения задачи и освобождает всю ранее выделенную память для хранения данных.
Реализация всех перечисленных функций может быть выполнена по аналогии с ранее рассмотренными примерами и предоставляется читателю в качестве самостоятельного упражнения.
2. Функция ParallelHyperQuickSort. Функция производит параллельную быструю сортировку согласно рассмотренному алгоритму.
Пример 9.1.
(html, txt)
Функция LocalDataSort выполняет сортировку блока данных на каждом процессоре, используя последовательный алгоритм быстрой сортировки.
Функция PivotDistribution определяет ведущий элемент и рассылает его значение всем процессорам.
Функция GetProcDataDivisionPos выполняет разделение блока данных относительно ведущего элемента. Ее результатом является целое число, обозначающее позицию элемента на границе двух блоков.
Функция DataMerge осуществляет слияние частей в один упорядоченный блок данных.
3. Функция PivotDistribution. Функция выбирает ведущий элемент и рассылает его все процессорам гиперкуба. Так как данные на процессорах отсортированы с самого начала, ведущий элемент выбирается как средний элемент блока данных.
Пример 9.2.
(html, txt)
Результаты вычислительных экспериментов
Эксперименты осуществлялись на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server.
Для оценки длительности ? базовой скалярной операции алгоритма сортировки проводилось решение задачи упорядочивания при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате выполненных экспериментов для величины ? было получено значение 10,41 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности
и пропускной способности ? соответственно 130 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, размер которого на данной платформе равен 8 байт (следовательно w=8).
Результаты вычислительных экспериментов приведены в табл. 9.2. Эксперименты выполнялись с использованием двух и четырех процессоров.
| 10000 | 0,001422 | 0,002210 | 0,643439 | 0,003270 | 0,434862 |
| 20000 | 0,002991 | 0,004428 | 0,675474 | 0,004596 | 0,650783 |
| 30000 | 0,004612 | 0,006745 | 0,683766 | 0,006873 | 0,671032 |
| 40000 | 0,006297 | 0,008033 | 0,783891 | 0,009107 | 0,691446 |
| 50000 | 0,008014 | 0,009770 | 0,820266 | 0,010840 | 0,739299 |
Рис. 9.1. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма пузырьковой сортировки
Как можно заметить из приведенных результатов вычислительных экспериментов, параллельный вариант алгоритма сортировки работает медленнее исходного последовательного метода пузырьковой сортировки, т.к. объем передаваемых данных между процессорами является достаточно большим и сопоставим с количеством выполняемых вычислительных операций (и этот дисбаланс объема вычислений и сложности операций передачи данных увеличивается с ростом числа процессоров).
Вычислительные эксперименты для оценки эффективности параллельного варианта сортировки Шелла осуществлялись при тех же условиях, что и ранее выполненные (см. п. 9.3.6).
Результаты вычислительных экспериментов приведены в табл. 9.4. Эксперименты проводились с использованием двух и четырех процессоров. Время указано в секундах.
| 10000 | 0,001422 | 0,002959 | 0,480568 | 0,007509 | 0,189373 |
| 20000 | 0,002991 | 0,004557 | 0,656353 | 0,009826 | 0,304396 |
| 30000 | 0,004612 | 0,006118 | 0,753841 | 0,012431 | 0,371008 |
| 40000 | 0,006297 | 0,008461 | 0,744238 | 0,017009 | 0,370216 |
| 50000 | 0,008014 | 0,009920 | 0,807863 | 0,019419 | 0,412689 |
Рис. 9.4. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма сортировки Шелла
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (9.7) приведено в таблице 9.5 и на рис. 9.5.
| 1000 | 0,002684 | 0,002959 | 0,002938 | 0,007509 |
| 20000 | 0,004872 | 0,004557 | 0,004729 | 0,009826 |
| 30000 | 0,007100 | 0,006118 | 0,006538 | 0,012431 |
| 40000 | 0,009353 | 0,008461 | 0,008361 | 0,017009 |
| 50000 | 0,011625 | 0,009920 | 0,010193 | 0,019419 |
Рис. 9.5. График зависимости экспериментального и теоретического времени проведения эксперимента на двух процессорах от объема исходных данных
Вычислительные эксперименты для оценки эффективности параллельного варианта быстрой сортировки производились при тех же условиях, что и ранее выполненные эксперименты (см. п. 9.3.6).
Результаты вычислительных экспериментов приведены в табл. 9.6. Эксперименты проводились с использованием двух и четырех процессоров. Время указано в секундах.
| 10000 | 0,001422 | 0,001521 | 0,934911 | 0,003434 | 0,414094 |
| 20000 | 0,002991 | 0,002234 | 1,338854 | 0,004094 | 0,730581 |
| 30000 | 0,004612 | 0,003080 | 1,497403 | 0,005088 | 0,906447 |
| 40000 | 0,006297 | 0,004363 | 1,443273 | 0,005906 | 1,066204 |
| 50000 | 0,008014 | 0,005486 | 1,460809 | 0,006635 | 1,207837 |
Как можно заметить по результатам вычислительных экспериментов, параллельный алгоритм быстрой сортировки уже позволяет получить ускорение при решении задачи упорядочивания данных.
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (9.11) приведено в таблице 9.7 и на рис. 9.8.
Рис. 9.7. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма быстрой сортировки
| 10000 | 0,001280 | 0,001521 | 0,001735 | 0,003434 |
| 20000 | 0,002265 | 0,002234 | 0,002321 | 0,004094 |
| 30000 | 0,003289 | 0,003080 | 0,002928 | 0,005088 |
| 40000 | 0,004338 | 0,004363 | 0,003547 | 0,005906 |
| 50000 | 0,005407 | 0,005486 | 0,004175 | 0,006635 |
Рис. 9.8. График зависимости экспериментального и теоретического времени проведения эксперимента на двух процессорах от объема исходных данных
Вычислительные эксперименты для оценки эффективности параллельного варианта обобщенной быстрой сортировки производились при тех же условиях, что и ранее выполненные (см. п. 9.3.6).
Результаты вычислительных экспериментов даны в табл. 9.8. Эксперименты проводились с использованием двух и четырех процессоров. Время указано в секундах.
| 10000 | 0,001422 | 0,001485 | 0,957576 | 0,002898 | 0,490683 |
| 20000 | 0,002991 | 0,002180 | 1,372018 | 0,003770 | 0,793369 |
| 30000 | 0,004612 | 0,003077 | 1,498863 | 0,004451 | 1,036172 |
| 40000 | 0,006297 | 0,003859 | 1,631770 | 0,004721 | 1,333828 |
| 50000 | 0,008014 | 0,005041 | 1,589764 | 0,005242 | 1,528806 |
Рис. 9.9. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма обобщенной быстрой сортировки
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (9.12) приведено в таблице 9.9 и на рис. 9.10.
| 10000 | 0,001281 | 0,001485 | 0,001735 | 0,002898 |
| 20000 | 0,002265 | 0,002180 | 0,002322 | 0,003770 |
| 30000 | 0,003289 | 0,003077 | 0,002928 | 0,004451 |
| 40000 | 0,004338 | 0,003859 | 0,003547 | 0,004721 |
| 50000 | 0,005407 | 0,005041 | 0,004176 | 0,005242 |
Рис. 9.10. График зависимости экспериментального и теоретического времени проведения эксперимента на четырех процессорах от объема исходных данных
Вычислительные эксперименты для оценки эффективности параллельного варианта сортировки с использованием регулярного набора образцов осуществлялись при тех же условиях, что и ранее выполненные (см. п. 9.3.6).
Результаты вычислительных экспериментов даны в табл. 9.10. Эксперименты проводились с использованием двух и четырех процессоров. Время указано в секундах.
| 10000 | 0,001422 | 0,001513 | 0,939855 | 0,001166 | 1,219554 |
| 20000 | 0,002991 | 0,002307 | 1,396489 | 0,002081 | 1,437290 |
| 30000 | 0,004612 | 0,003168 | 1,455808 | 0,003099 | 1,488222 |
| 40000 | 0,006297 | 0,004542 | 1,386394 | 0,003819 | 1,648861 |
| 50000 | 0,008014 | 0,005503 | 1,456297 | 0,004370 | 1,833867 |
Рис. 9.12. Зависимость ускорения от количества процессоров при выполнении параллельного алгоритма сортировки с использованием регулярного набора образцов
| 10000 | 0,001533 | 0,001513 | 0,001762 | 0,001166 |
| 20000 | 0,002569 | 0,002307 | 0,002375 | 0,002081 |
| 30000 | 0,003645 | 0,003168 | 0,003007 | 0,003099 |
| 40000 | 0,004747 | 0,004542 | 0,003652 | 0,003819 |
| 50000 | 0,005867 | 0,005503 | 0,004307 | 0,004370 |
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (9.17) приведено в таблице 9.11 и на рис. 9.13.
Рис. 9.13. График зависимости экспериментального и теоретического времени проведения эксперимента на четырех процессорах от объема исходных данных
Выполните реализацию параллельного алгоритма пузырьковой
Выполните реализацию параллельного алгоритма пузырьковой сортировки. Проведите эксперименты. Постройте теоретические оценки с учетом тех операций пересылок данных, которые использовались при реализации, и параметров вычислительной системы. Сравните получаемые теоретические оценки с результатами экспериментов.Выполните реализацию параллельного алгоритма быстрой сортировки по одной из приведенных схем. Определите значения параметров латентности, пропускной способности и времени выполнения базовой операции для используемой вычислительной системы и получите оценки показателей ускорения и эффективности для реализованного метода параллельных вычислений.Разработайте параллельную схему вычислений для широко известного алгоритма сортировки слиянием (подробное описание метода может быть получено, например, в работах [[26]] или [[50]]). Выполните реализацию разработанного алгоритма и постройте все необходимые теоретические оценки сложности метода.
Алгоритм Кернигана – Лина
В алгоритме Кернигана – Лина (the Kernighan – Lin algorithm) используется несколько иной подход для решения проблемы оптимального разбиения графа – предполагается, что некоторое начальное разбиение графа уже существует, затем имеющееся приближение улучшается в течение некоторого количества итераций. Применяемый способ улучшения в алгоритме Кернигана – Лина состоит в обмене вершинами между подмножествами имеющегося разбиения графа (см. рис. 10.16). Для формирования требуемого количества частей графа может быть использована, как и ранее, рекурсивная процедура деления пополам.
Общая схема одной итерации алгоритма Кернигана – Лина может быть представлена следующим образом.
Алгоритм 10.3. Общая схема алгоритма Кернигана – Лина
Формирование множества пар вершин для перестановки. Из вершин, которые еще не были переставлены на данной итерации, формируются все возможные пары (в парах должно присутствовать по одной вершине из каждой части имеющегося разбиения графа).Построение новых вариантов разбиения графа. Каждая пара, подготовленная на шаге 1, поочередно используется для обмена вершин между частями имеющегося разбиения графа для получения множества новых вариантов деления.Выбор лучшего варианта разбиения графа. Для сформированного на шаге 2 множества новых делений графа выбирается лучший вариант. Этот вариант далее фиксируется как новое текущее разбиение графа, а соответствующая выбранному варианту пара вершин отмечается как использованная на текущей итерации алгоритма.Проверка использования всех вершин. При наличии в графе вершин, еще не использованных при перестановках, выполнение итерации алгоритма снова продолжается с шага 1. Если же перебор вершин графа завершен, далее следует шаг 5.Выбор наилучшего варианта разбиения графа. Среди всех разбиений графа, полученных на шаге 3 проведенных итераций, выбирается (и фиксируется) наилучший вариант разбиения графа.
Поясним дополнительно, что на шаге 2 итерации алгоритма перестановка вершин каждой очередной пары осуществляется для одного и того же разбиения графа, выбранного до начала выполнения итерации или определенного на шаге 3. Общее количество выполняемых итераций, как правило, фиксируется заранее и является параметром алгоритма (за исключением случая остановки при отсутствии улучшения разбиения на очередной итерации).
Рис. 10.16. Пример перестановки двух вершин (выделены серым) в методе Кернигана – Лина
Анализ эффективности параллельных вычислений
Выполним анализ эффективности параллельного алгоритма Флойда, обеспечивающего поиск всех кратчайших путей. Как и ранее, проведем этот анализ в два этапа. На первом оценим порядок вычислительной сложности алгоритма, затем на втором этапе уточним полученные оценки и учтем трудоемкость выполнения коммуникационных операций.
Общая трудоемкость последовательного алгоритма, как уже отмечалось ранее, имеет порядок сложности n3. Для параллельного алгоритма на отдельной итерации каждый процессор выполняет обновление элементов матрицы А. Всего в подзадачах n2/p таких элементов, число итераций алгоритма равно n – таким образом, показатели ускорения и эффективности параллельного алгоритма Флойда имеют вид:
|
(10.1) |
Следовательно, общий анализ сложности дает идеальные показатели эффективности параллельных вычислений. Для уточнения полученных соотношений введем в полученные выражения время выполнения базовой операции выбора минимального значения и учтем затраты на выполнение операций передачи данных между процессорами.
Коммуникационная операция, выполняемая на каждой итерации алгоритма Флойда, состоит в передаче одной из строк матрицы А всем процессорам вычислительной системы. Как уже показывалось ранее, такая операция может быть выполнена за ?log2p? шагов. С учетом количества итераций алгоритма Флойда при использовании модели Хокни общая длительность выполнения коммуникационных операций может быть определена при помощи следующего выражения
|
(10.2) |
где, как и ранее, – латентность сети передачи данных, ? – пропускная способность сети, а w есть размер элемента матрицы в байтах.
С учетом полученных соотношений общее время выполнения параллельного алгоритма Флойда может быть определено следующим образом:
|
(10.3) |
где ? есть время выполнения операции выбора минимального значения.
Общий анализ сложности параллельного алгоритма Прима для нахождения минимального охватывающего дерева дает идеальные показатели эффективности параллельных вычислений:
|
(10.4) |
При этом следует отметить, что в ходе параллельных вычислений идеальная балансировка вычислительной нагрузки процессоров может быть нарушена. В зависимости от вида исходного графа G количество выбранных вершин в охватывающем дереве на разных процессорах может оказаться различным, и распределение вычислений между процессорами станет неравномерным (вплоть до отсутствия вычислительной нагрузки на отдельных процессорах). Однако такие предельные ситуации нарушения балансировки в общем случае возникают достаточно редко, а организация динамического перераспределения вычислительной нагрузки между процессорами в ходе вычислений является интересной, но одновременно и очень сложной задачей.
Для уточнения полученных показателей эффективности параллельных вычислений оценим более точно количество вычислительных операций алгоритма и учтем затраты на выполнение операций передачи данных между процессорами.
При выполнении вычислений на отдельной итерации параллельного алгоритма Прима каждый процессор определяет номер ближайшей вершины из Vj до охватывающего дерева и осуществляет корректировку расстояний di после расширения МОД. Количество выполняемых операций в каждой из этих вычислительных процедур ограничивается сверху числом вершин, имеющихся на процессорах, т.е. величиной ?n/p?. Как результат, с учетом общего количества итераций n время выполнения вычислительных операций параллельного алгоритма Прима может быть оценено при помощи соотношения:
|
(10.5) |
(здесь, как и ранее, ? есть время выполнения одной элементарной скалярной операции).
Операция сбора данных от всех процессоров на одном из процессоров может быть произведена за ?log2p? итераций, при этом общая оценка длительности выполнения передачи данных определяется выражением (более подробное рассмотрение данной коммуникационной операции содержится в лекции 3):
|
(10.6) |
где – латентность сети передачи данных, ? – пропускная способность сети, а w есть размер одного пересылаемого элемента данных в байтах (коэффициент 3 в выражении соответствует числу передаваемых значений между процессорами – длина минимальной дуги и номера двух вершин, которые соединяются этой дугой).
Коммуникационная операция передачи данных от одного процессора всем процессорам вычислительной системы также может быть выполнена за ?log2p? итераций при общей оценке времени выполнения вида:
|
(10.7) |
С учетом всех полученных соотношений общее время выполнения параллельного алгоритма Прима составляет:
|
(10.8) |
Деление с учетом связности
С самых общих позиций понятно, что при разделении графа информационная зависимость между разделенными подграфами будет меньше, если соседние вершины (вершины, между которыми имеются дуги) будут находиться в одном подграфе. Алгоритм деления графов с учетом связности (the levelized nested dissection algorithm) пытается достичь этого, последовательно добавляя к формируемому подграфу соседей. На каждой итерации алгоритма происходит разделение графа на 2 части. Таким образом, разделение графа на требуемое число частей достигается путем рекурсивного применения алгоритма.
Общая схема алгоритма может быть описана при помощи следующего набора правил.
Алгоритм 10.2. Общая схема выполнения алгоритма деления графов с учетом связности
Iteration = 0.Присвоение номера Iteration произвольной вершине графа.Присвоение ненумерованным соседям вершин с номером Iteration номера Iteration + 1.Iteration = Iteration + 1.Если еще есть неперенумерованные соседи, то переход на шаг 3.Разделение графа на 2 части в порядке нумерации.
Для минимизации информационных зависимостей имеет смысл в качестве начальной выбирать граничную вершину. Поиск такой вершины можно осуществить методом, близким к рассмотренной схеме. Так, перенумеровав вершины графа в соответствии с алгоритмом 10.2 (начиная нумерацию из произвольной вершины), мы можем взять любую вершину с максимальным номером. Как нетрудно убедиться, она будет граничной.
Пример работы алгоритма приведен на рис. 10.15. Цифрами показаны номера, которые получили вершины в процессе разделения. Сплошной линией показана граница, разделяющая 2 подграфа. Также на рисунке показано лучшее решение (пунктирная линия). Очевидно, что полученное алгоритмом разбиение далеко от оптимального, так как в приведенном примере есть решение только с тремя пересеченными ребрами вместо пяти.
Рис. 10.15. Пример работы алгоритма деления графов с учетом связности
Деление сети с использованием кривых Пеано
Одним из недостатков предыдущих графических методов является то, что при каждой бисекции эти методы учитывают только одну размерность. Таким образом, схемы, учитывающие больше размерностей, могут обеспечить лучшее разбиение.
Один из таких методов упорядочивает элементы в соответствии с позициями центров их масс вдоль кривых Пеано. Кривые Пеано – это кривые, полностью заполняющие области больших размерностей (например, квадрат или куб). Применение таких кривых обеспечивает близость точек фигуры, которые соответствуют точкам, близким на кривой. После получения списка элементов сети, упорядоченного в зависимости от расположения на кривой, достаточно разделить список на необходимое число частей в соответствии с установленным порядком. Получаемый в результате такого подхода метод носит в литературе наименование алгоритма деления сети с использованием кривых Пеано (the space-filling curve technique). Подробнее о методе можно прочитать в работах [[56], [58], [61]].
Рис. 10.14. Пример разделения сети на 3 части с использованием кривых Пеано
Геометрические методы
Геометрические методы (см., например, [[21], [35], [44], [53], [55], [58], [61], [65]]) выполняют разбиение сетей, основываясь исключительно на координатной информации об узлах сети. Так как эти методы не принимают во внимание информацию о связности элементов сети, они не могут явно привести к минимизации суммарного веса граничных ребер (в терминах графа, соответствующего сети). Для минимизации межпроцессорных коммуникаций геометрические методы оптимизируют некоторые вспомогательные показатели (например, длину границы между разделенными участками сети).
Обычно геометрические методы не требуют большого объема вычислений, однако качество их разбиения уступает методам, принимающим во внимание связность элементов сети.
Комбинаторные методы
В отличие от геометрических методов, комбинаторные алгоритмы (см., например, [[36], [67]]) обычно оперируют не с сетью, а с графом, построенным для этой сети. Соответственно, в отличие от геометрических схем, комбинаторные методы не принимают во внимание информацию о близости расположения элементов сети друг относительно друга, руководствуясь только смежностью вершин графа. Комбинаторные методы обычно обеспечивают более сбалансированное разбиение и меньшее информационное взаимодействие полученных подсетей. Однако комбинаторные методы имеют тенденцию работать существенно дольше, чем их геометрические аналоги.
Контрольные вопросы
Приведите определение графа. Какие основные способы используются для задания графов?В чем состоит задача поиска всех кратчайших путей?Приведите общую схему алгоритма Флойда. Какова трудоемкость алгоритма?В чем состоит способ распараллеливания алгоритма Флойда?В чем заключается задача нахождения минимального охватывающего дерева? Приведите пример использования задачи на практике.Приведите общую схему алгоритма Прима. Какова трудоемкость алгоритма?В чем состоит способ распараллеливания алгоритма Прима?В чем отличие геометрических и комбинаторных методов разделения графа? Какие методы являются более предпочтительными? Почему?Приведите описание метода покоординатного разбиения и алгоритма разделения с учетом связности. Какой из этих методов является более простым для реализации?
Краткий обзор лекции
В лекции рассмотрен ряд алгоритмов для решения типовых задач обработки графов. Кроме того, приведен обзор методов разделения графа.
В подразделе 10.1 представлен алгоритм Флойда (the Floyd algorithm) – дается общая вычислительная схема последовательного варианта метода, обсуждаются способы его распараллеливания, проводится анализ эффективности получаемых параллельных вычислений, рассматривается программная реализация метода и приводятся результаты вычислительных экспериментов. Используемый подход к распараллеливанию алгоритма Флойда состоит в разделении вершин графа между процессорами, а необходимое при этом информационное взаимодействие состоит в передаче одной строки матрицы смежности от одного процессора всем процессорам вычислительной системы на каждой итерации метода.
В подразделе 10.2 рассматривается алгоритм Прима (the Prim algorithm) для решения задачи поиска минимального охватывающего дерева (остова) неориентированного взвешенного графа. Остовом графа называют связный подграф без циклов (дерево), который содержит все вершины исходного графа и ребра, имеющие минимальный суммарный вес. Для алгоритма дается общее описание его исходного последовательного варианта, определяются возможные способы его параллельного выполнения, теоретические оценки ускорения и эффективности параллельных вычислений; также рассматриваются результаты проведенных вычислительных экспериментов. Параллельный вариант алгоритма Прима, как и в предыдущем случае, основывается на разделении вершин графа между процессорами при несколько большем объеме информационных взаимодействий – на каждой итерации алгоритма необходимой является операция сбора данных на одном процессоре и последующая рассылка номера выбранной вершины графа всем процессорам вычислительной системы.
Рассматриваемая в подразделе 10.3 задача оптимального разделения графов является важной для многих научных исследований, использующих параллельные вычисления. Для примера в подразделе приведен общий способ перехода от двумерной или трехмерной сети, моделирующей процесс вычислений, к соответствующему ей графу. Для решения задачи разбиения графов были рассмотрены геометрические методы, использующие при разделении сетей только координатную информацию об узлах сети, и комбинаторные алгоритмы, руководствующиеся смежностью вершин графа. К числу рассмотренных геометрических методов относятся покоординатное разбиение (the coordinate nested dissection method), рекурсивный инерционный метод деления пополам (the recursive inertial bisection method), деление сети с использованием кривых Пеано (the space-filling curve techniques). К числу рассмотренных комбинаторных алгоритмов относятся деление с учетом связности (the levelized nested dissection) и алгоритм Кернигана – Лина (the Kernighan – Lin algorithm). Для сопоставления рассмотренных подходов приводится общая сравнительная характеристика алгоритмов по времени выполнения, точности получаемого решения, возможностям для распараллеливания и т.п.
Масштабирование и распределение подзадач по процессорам
Как правило, число доступных процессоров p существенно меньше, чем число базовых задач n2 (p<<n2). Возможный способ укрупнения вычислений состоит в использовании ленточной схемы разбиения матрицы A – такой подход соответствует объединению в рамках одной базовой подзадачи вычислений, связанных с обновлением элементов одной или нескольких строк (горизонтальное разбиение) или столбцов (вертикальное разбиение) матрицы A. Эти два типа разбиения практически равноправны – учитывая дополнительный момент, что для алгоритмического языка C массивы располагаются по строкам, будем рассматривать далее только разбиение матрицы A на горизонтальные полосы.
Следует отметить, что при таком способе разбиения данных на каждой итерации алгоритма Флойда потребуется передавать между подзадачами только элементы одной из строк матрицы A. Для оптимального выполнения подобной коммуникационной операции топология сети должна обеспечивать эффективное представление структуры сети передачи данных в виде гиперкуба или полного графа.
По определению количество базовых подзадач всегда соответствует числу имеющихся процессоров, и, тем самым, проблема масштабирования для параллельного алгоритма не возникает.
Распределение подзадач между процессорами должно учитывать характер выполняемых в алгоритме Прима коммуникационных операций. Для оптимальной реализации требуемых информационных взаимодействий между базовыми подзадачами топология сети передачи данных должна обеспечивать эффективное представление в виде гиперкуба или полного графа.
Метод рекурсивного деления пополам
Для решения задачи разбиения графа можно рекурсивно применить метод бинарного деления (the binary bisection method), при котором на первой итерации граф разделяется на две равные части, далее на втором шаге каждая из полученных частей также разбивается на две части и т. д. В данном подходе для разбиения графа на k частей необходимо log2k уровней рекурсии и выполнение k-1 деления пополам. В случае когда требуемое количество разбиений k не является степенью двойки, каждое деление пополам необходимо осуществлять в соответствующем соотношении.
Поясним схему работы метода деления пополам на примере разделения графа на рис. 10.11 на 5 частей. Сначала граф следует разделить на 2 части в оотношении 2:3 (непрерывная линия), затем правую часть разбиения – в отношении 1:3 (штриховая линия), после этого осталось разделить 2 крайние подобласти слева и справа в отношении 1:1 (штрих-пунктир).
Рис. 10.11. Пример разбиения графа на 5 частей методом рекурсивного деления пополам
Обзор литературы
Дополнительная информация по алгоритмам Флойда и Прима может быть получена, например, в [[26]].
Подробное рассмотрение вопросов, связанных с проблемой разделения графов, содержится в работах [[21], [36], [37], [44], [53], [55], [58], [61], [65], [67]].
Параллельные алгоритмы разделения графов рассматриваются в [[20], [38], [44], [48], [49], [65], [74]].
Параллельные методы на графах
Математические модели в виде графов широко используются при моделировании разнообразных явлений, процессов и систем. Как результат, многие теоретические и реальные прикладные задачи могут быть решены при помощи тех или иных процедур анализа графовых моделей. Среди множества этих процедур может быть выделен некоторый определенный набор типовых алгоритмов обработки графов. Рассмотрению вопросов теории графов, алгоритмов моделирования, анализу и решению задач на графах посвящено достаточно много различных изданий, в качестве возможного руководства по данной тематике может быть рекомендована работа [[26]].
Пусть G есть граф
G=(V,R),
для которого набор вершин Vi, 0in, задается множеством V, а список дуг графа
определяется множеством R. В общем случае дугам графа могут приписываться некоторые числовые характеристики (веса) wj, 0jm (взвешенный граф). Пример взвешенного графа приведен на рис. 10.1.
Рис. 10.1. Пример взвешенного ориентированного графа
Известны различные способы задания графов. При малом количестве дуг в графе (т. е. m<<n2) целесообразно использовать для определения графов списки, перечисляющие имеющиеся в графах дуги. Представление достаточно плотных графов, для которых почти все вершины соединены между собой дугами (т. е. m~n2), может быть эффективно обеспечено при помощи матрицы смежности:
A=(aij), 1i, jn,
ненулевые значения элементов которой соответствуют дугам графа:
(для обозначения отсутствия ребра между вершинами в матрице смежности на соответствующей позиции используется знак бесконечности, при вычислениях знак бесконечности может быть заменен, например, на любое отрицательное число). Так, например, матрица смежности, соответствующая графу на рис. 10.1, приведена на рис. 10.2.
Рис. 10.2. Матрица смежности для графа с рис. 10.1
Как положительный момент такого способа представления графов можно отметить, что использование матрицы смежности позволяет применять при реализации вычислительных процедур анализа графов матричные алгоритмы обработки данных.
Далее мы рассмотрим способы параллельной реализации алгоритмов на графах на примере задачи поиска кратчайших путей между всеми парами пунктов назначения и задачи выделения минимального охватывающего дерева (остова) графа. Кроме того, мы рассмотрим задачу оптимального разделения графов, широко используемую для организации параллельных вычислений. Для представления графов при рассмотрении всех перечисленных задач будут применяться матрицы смежности.
Покоординатное разбиение
Покоординатное разбиение (the coordinate nested dissection) – это метод, основанный на рекурсивном делении пополам сети по наиболее длинной стороне. В качестве иллюстрации на рис. 10.12
показан пример сети, при разделении которой именно такой способ разбиения дает существенно меньшее количество информационных связей между разделенными частями, по сравнению со случаем, когда сеть делится по меньшей (вертикальной) стороне.
Рис. 10.12. Пример разделения сети графическим методом по наибольшей размерности (граница раздела показана жирной линией)
Общая схема выполнения метода состоит в следующем. Сначала вычисляются центры масс элементов сети. Полученные точки проектируются на ось, соответствующую наибольшей стороне разделяемой сети. Таким образом мы получаем упорядоченный список всех элементов сети. Делением списка пополам (возможно, в нужной пропорции) мы получаем требуемую бисекцию. Аналогичным способом полученные фрагменты разбиения рекурсивно делятся на нужное число частей.
Метод координатного вложенного разбиения работает очень быстро и требует небольшого количества оперативной памяти. Однако получаемое разбиение уступает по качеству более сложным и вычислительно трудоемким методам. Кроме того, в случае сложной структуры сети алгоритм может получать разбиение с несвязанными подсетями.
Последовательный алгоритм Флойда
Для поиска минимальных расстояний между всеми парами пунктов назначения Флойд предложил алгоритм, сложность которого имеет порядок n3. В общем виде данный алгоритм может быть представлен следующим образом:
Алгоритм 10.1. Общая схема алгоритма Флойда
// Алгоритм 10.1 // Последовательный алгоритм Флойда for (k = 0; k < n; k++) for (i = 0; i < n; i++) for (j = 0; j < n; j++) A[i, j] = min(A[i, j], A[i, k] + A[k, j]);
(реализация операции выбора минимального значения min должна учитывать способ указания в матрице смежности несуществующих дуг графа). Как можно заметить, в ходе выполнения алгоритма матрица смежности A изменяется, после завершения вычислений в матрице A будет храниться требуемый результат – длины минимальных путей для каждой пары вершин исходного графа.
Дополнительная информация и доказательство правильности алгоритма Флойда могут быть получены, например, в работе [[26]].
Последовательный алгоритм Прима
Алгоритм начинает работу с произвольной вершины графа, выбираемой в качестве корня дерева, и в ходе последовательно выполняемых итераций расширяет конструируемое дерево до МОД. Пусть VT есть множество вершин, уже включенных алгоритмом в МОД, а величины di, 1in, характеризуют дуги минимальной длины от вершин, еще не включенных в дерево, до множества VT, т.е.
(если для какой-либо вершины iVT не существует ни одной дуги в VT, значение di устанавливается равным ?). В начале работы алгоритма выбирается корневая вершина МОД s и полагается VT={s}, ds=0.
Действия, выполняемые на каждой итерации алгоритма Прима, состоят в следующем:
определяются значения величин di для всех вершин, еще не включенных в состав МОД;выбирается вершина t графа G, имеющая дугу минимального веса до множества VTt:dt, iVT;вершина t включается в VT.
После выполнения n-1 итерации метода МОД будет сформировано. Вес этого дерева может быть получен при помощи выражения
Трудоемкость нахождения МОД характеризуется квадратичной зависимостью от числа вершин графа T1~n2.
Постановка задачи оптимального разделения графов
Пусть дан взвешенный неориентированный граф G=(V,E), каждой вершине ?V и каждому ребру eE которого приписан вес. Задача оптимального разделения графа состоит в разбиении его вершин на непересекающиеся подмножества с максимально близкими суммарными весами вершин и минимальным суммарным весом ребер, проходящих между полученными подмножествами вершин.
Следует отметить возможную противоречивость указанных критериев разбиения графа – равновесность подмножеств вершин может не соответствовать минимальности весов граничных ребер и наоборот. В большинстве случаев необходимым является выбор того или иного компромиссного решения. Так, в случае невысокой доли коммуникаций может оказаться эффективной оптимизация веса ребер только среди решений, обеспечивающих оптимальное разбиение множества вершин по весу.
Далее для простоты изложения учебного материала будем полагать веса вершин и ребер графа равными единице.
int A, int B)
// Программа 10.1. — Алгоритм Флойда
int ProcRank; // Ранг текущего процесса
int ProcNum; // Количество процессов
// Функция вычисления минимума
int Min( int A, int B) {
int Result = (A < B) ? A : B;
if((A < 0) && (B >= 0)) Result = B;
if((B < 0) && (A >= 0)) Result = A;
if((A < 0) && (B < 0)) Result = -1;
return Result;
}
// Главная функция программы
int main(int argc, char* argv[]) {
int *pMatrix; // Матрица смежности
int Size; // Размер матрицы смежности
int *pProcRows; // Строки матрицы смежности текущего процесса
int RowNum; // Число строк для текущего процесса
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);
MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank);
// Инициализация данных
ProcessInitialization(pMatrix, pProcRows, Size, RowNum);
// Распределение данных между процессами
DataDistribution(pMatrix, pProcRows, Size, RowNum);
// Параллельный алгоритм Флойда
ParallelFloyd(pProcRows, Size, RowNum);
// Сбор результатов работы алгоритма
ResultCollection(pMatrix, pProcRows, Size, RowNum);
// Завершение вычислений процесса
ProcessTermination(pMatrix, pProcRows);
MPI_Finalize();
return 0;
}
|
// Программа 10.1. — Алгоритм Флойда int ProcRank; // Ранг текущего процесса int ProcNum; // Количество процессов // Функция вычисления минимума int Min(int A, int B) { int Result = (A < B) ? A : B; if((A < 0) && (B >= 0)) Result = B; if((B < 0) && (A >= 0)) Result = A; if((A < 0) && (B < 0)) Result = -1; return Result; } // Главная функция программы int main(int argc, char* argv[]) { int *pMatrix; // Матрица смежности int Size; // Размер матрицы смежности int *pProcRows; // Строки матрицы смежности текущего процесса int RowNum; // Число строк для текущего процесса MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &ProcNum); MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank); // Инициализация данных ProcessInitialization(pMatrix, pProcRows, Size, RowNum); // Распределение данных между процессами DataDistribution(pMatrix, pProcRows, Size, RowNum); // Параллельный алгоритм Флойда ParallelFloyd(pProcRows, Size, RowNum); // Сбор результатов работы алгоритма ResultCollection(pMatrix, pProcRows, Size, RowNum); // Завершение вычислений процесса ProcessTermination(pMatrix, pProcRows); MPI_Finalize(); return 0; } |
| Пример 10.1. |
| Закрыть окно |
// Функция для рассылки строки всем процессам
void RowDistribution(int *pProcRows, int Size, int RowNum, int k,
int *pRow) {
int ProcRowRank; // Ранг процесса, которому принадлежит k-я строка
int ProcRowNum; // Номер k-й строки в полосе матрицы
// Нахождение ранга процесса – владельца k-й строки
int RestRows = Size;
int Ind = 0;
int Num = Size / ProcNum;
for(ProcRowRank = 1; ProcRowRank < ProcNum + 1; ProcRowRank ++) {
if(k < Ind + Num ) break;
RestRows -= Num;
Ind += Num;
Num = RestRows / (ProcNum - ProcRowRank);
}
ProcRowRank = ProcRowRank - 1;
ProcRowNum = k - Ind;
if(ProcRowRank == ProcRank)
// Копирование строки в массив pRow
copy(&pProcRows[ProcRowNum * Size], &pProcRows[(ProcRowNum + 1) *
Size], pRow);
// Распределение k-й строки между процессами
MPI_Bcast(pRow, Size, MPI_INT, ProcRowRank, MPI_COMM_WORLD);
}
|
// Функция для рассылки строки всем процессам void RowDistribution(int *pProcRows, int Size, int RowNum, int k, int *pRow) { int ProcRowRank; // Ранг процесса, которому принадлежит k-я строка int ProcRowNum; // Номер k-й строки в полосе матрицы // Нахождение ранга процесса – владельца k-й строки int RestRows = Size; int Ind = 0; int Num = Size / ProcNum; for(ProcRowRank = 1; ProcRowRank < ProcNum + 1; ProcRowRank ++) { if(k < Ind + Num ) break; RestRows -= Num; Ind += Num; Num = RestRows / (ProcNum - ProcRowRank); } ProcRowRank = ProcRowRank - 1; ProcRowNum = k - Ind; if(ProcRowRank == ProcRank) // Копирование строки в массив pRow copy(&pProcRows[ProcRowNum * Size], &pProcRows[(ProcRowNum + 1) * Size], pRow); // Распределение k-й строки между процессами MPI_Bcast(pRow, Size, MPI_INT, ProcRowRank, MPI_COMM_WORLD); } |
| Пример 10.2. |
| Закрыть окно |
Программная реализация
Представим возможный вариант параллельной реализации алгоритма Флойда. При этом описание отдельных модулей не приводится, если их отсутствие не оказывает влияния на понимание общей схемы параллельных вычислений.
1. Главная функция программы. Реализует логику работы алгоритма, последовательно вызывает необходимые подпрограммы.
Пример 10.1.
(html, txt)
Функция Min вычисляет меньшее из двух целых чисел, учитывая применяемый способ задания несуществующих дуг в матрице смежности (в рассматриваемой реализации для этого используется значение -1).
Функция ProcessInitialization определяет исходные данные решаемой задачи (количество вершин графа), выделяет память для хранения данных, осуществляет ввод матрицы смежности (или формирует ее при помощи какого-либо датчика случайных числе).
Функция DataDistribution распределяет исходные данные между процессами. Каждый процесс получает горизонтальную полосу матрицы смежности.
Функция ResultCollection осуществляет сбор со всех процессов горизонтальных полос результирующей матрицы кратчайших расстояний между любыми парами вершин графа.
Функция ProcessTermination выполняет необходимый вывод результатов решения задачи и освобождает всю ранее выделенную память для хранения данных.
Реализация всех перечисленных функций может быть выполнена по аналогии с ранее рассмотренными примерами и предоставляется читателю в качестве самостоятельного упражнения.
2. Функция ParallelFloyd. Данная функция осуществляет итеративное изменение матрицы смежности в соответствии с алгоритмом Флойда.
// Параллельный алгоритм Флойда void ParallelFloyd(int *pProcRows, int Size, int RowNum) { int *pRow = new int[Size]; int t1, t2;
for(int k = 0; k < Size; k++) { // Распределение k-й строки среди процессов RowDistribution(pProcRows, Size, RowNum, k, pRow);
// Обновление элементов матрицы смежности for(int i = 0; i < RowNum; i++) for(int j = 0; j < Size; j++) if( (pProcRows[i * Size + k] != -1) && (pRow [j]!= -1)){ t1 = pProcRows[i * Size + j]; t2 = pProcRows[i * Size + k] + pRow[j]; pProcRows[i * Size + j] = Min(t1, t2); } }
delete []pRow; }
3. Функция RowDistribution. Данная функция рассылает k-ю строку матрицы смежности всем процессам программы.
Пример 10.2.
(html, txt)
Разделение вычислений на независимые части
Оценим возможности параллельного выполнения рассмотренного алгоритма нахождения минимально охватывающего дерева.
Итерации метода должны выполняться последовательно и, тем самым, не могут быть распараллелены. С другой стороны, выполняемые на каждой итерации алгоритма действия являются независимыми и могут реализовываться одновременно. Так, например, определение величин di может осуществляться для каждой вершины графа в отдельности, нахождение дуги минимального веса может быть реализовано по каскадной схеме и т.д.
Распределение данных между процессорами вычислительной системы должно обеспечивать независимость перечисленных операций алгоритма Прима. В частности, это может быть реализовано, если каждая вершина графа располагается на процессоре вместе со всей связанной с вершиной информацией. Соблюдение данного принципа приводит к тому, что при равномерной загрузке каждый процессор pj, 1jp, должен содержать:
набор вершин
соответствующий этому набору блок из k величин
вертикальную полосу матрицы смежности графа G из k соседних столбцов
общую часть набора Vj и формируемого в процессе вычислений множества вершин VT.
Как итог можем заключить, что базовой подзадачей в параллельном алгоритме Прима может служить процедура вычисления блока значений ?j для вершин Vj матрицы смежности A графа G.
Как следует из общей схемы алгоритма Флойда, основная вычислительная нагрузка при решении задачи поиска кратчайших путей состоит в выполнении операции выбора минимальных значений (см. Алгоритм 10.1). Данная операция является достаточно простой, и ее распараллеливание не приведет к заметному ускорению вычислений. Более эффективный способ организации параллельных вычислений может состоять в одновременном выполнении нескольких операций обновления значений матрицы A.
Покажем корректность такого способа организации параллелизма. Для этого нужно доказать, что операции обновления значений матрицы A на одной и той же итерации внешнего цикла k могут выполняться независимо. Иными словами, следует показать, что на итерации k не происходит изменения элементов Aik и Akj ни для одной пары индексов (i, j). Рассмотрим выражение, по которому происходит изменение элементов матрицы A:
Aij min (Aij, Aik + Akj).
Для i=k получим
Akj min (Akj, Akk + Akj),
но тогда значение Akj не изменится, т.к. Akk=0.
Для j=k выражение преобразуется к виду
Aik min (Aik, Aik + Akk),
что также показывает неизменность значений Aik. Как результат, необходимые условия для организации параллельных вычислений обеспечены, и, тем самым, в качестве базовой подзадачи может быть использована операция обновления элементов матрицы A (для указания подзадач будем применять индексы обновляемых в подзадачах элементов).
Рекурсивный инерционный метод деления пополам
Предыдущая схема могла производить разбиение сети только по линии, перпендикулярной одной из координатных осей. Во многих случаях такое ограничение оказывается критичным для построения качественного разбиения. Достаточно повернуть сеть на рис. 10.12
под острым углом к координатным осям (см. рис. 10.13), чтобы убедиться в этом. Для минимизации границы между подсетями желательна возможность проведения линии разделения с любым требуемым углом поворота. Возможный способ определения угла поворота, используемый в рекурсивном инерционном методе деления пополам (the recursive inertial bisection), состоит в использовании главной инерционной оси (см., например, [[62]]), считая элементы сети точечными массами. Линия бисекции, ортогональная полученной оси, как правило, дает границу наименьшей длины.
Рис. 10.13. Пример разделения сети методом рекурсивной инерционной бисекции. Стрелкой показана главная инерционная ось
Результаты вычислительных экспериментов
Эксперименты проводились на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server.
Для оценки длительности ? базовой скалярной операции выбора минимального значения проводилось решение задачи поиска всех кратчайших путей при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате для величины ? было получено значение 7,14 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности и пропускной способности ? соответственно 130 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа int, размер которого на данной платформе равен 4 байта (следовательно, w=4).
Результаты вычислительных экспериментов приведены в таблице 10.1. Эксперименты выполнялись с использованием двух, четырех и восьми процессоров. Время указано в секундах.
| 1000 | 8,037 | 4,152 | 1,936 | 2,067 | 3,888 | 0,941 | 8,544 |
| 2000 | 59,812 | 30,323 | 1,972 | 15,375 | 3,890 | 8,058 | 7,423 |
| 3000 | 197,111 | 99,264 | 1,986 | 50,232 | 3,924 | 25,643 | 7,687 |
| 4000 | 461,785 | 232,507 | 1,986 | 117,220 | 3,939 | 69,946 | 6,602 |
| 5000 | 884,622 | 443,747 | 1,994 | 224,441 | 3,941 | 128,078 | 6,907 |
Рис. 10.4. Графики зависимости ускорения параллельного алгоритма Флойда от числа используемых процессоров при различном количестве вершин графа
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (10.3) приведено в таблице 10.2 и на рис. 10.5.
Рис. 10.5. Графики экспериментально установленного времени работы параллельного алгоритма Флойда и теоретической оценки в зависимости от количества вершин графа при использовании двух процессоров
Вычислительные эксперименты для оценки эффективности параллельного алгоритма Прима осуществлялись при тех же условиях, что и ранее выполненные (см. п. 10.1.7).
Для оценки длительности ? базовой скалярной операции проводилось решение задачи нахождения минимального охватывающего дерева при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины ? было получено значение 4,76 нсек. Все вычисления производились над числовыми значениями типа int, размер которого на данной платформе равен 4 байта (следовательно, w=4).
Результаты вычислительных экспериментов даны в таблице 10.3. Эксперименты проводились с использованием двух, четырех и восьми процессоров. Время указано в секундах.
| 1000 | 0,044 | 0,248 | 0,176 | 0,932 | 0,047 | 1,574 | 0,028 |
| 2000 | 0,208 | 0,684 | 0,304 | 1,800 | 0,115 | 2,159 | 0,096 |
| 3000 | 0,485 | 1,403 | 0,346 | 2,214 | 0,219 | 3,195 | 0,152 |
| 4000 | 0,873 | 1,946 | 0,622 | 3,324 | 0,263 | 5,431 | 0,161 |
| 5000 | 1,432 | 2,665 | 0,736 | 2,933 | 0,488 | 4,119 | 0,348 |
| 6000 | 2,189 | 2,900 | 0,821 | 4,291 | 0,510 | 7,737 | 0,283 |
| 7000 | 3,042 | 3,236 | 0,940 | 6,327 | 0,481 | 8,825 | 0,345 |
| 8000 | 4,150 | 4,462 | 0,930 | 6,993 | 0,593 | 10,390 | 0,399 |
| 9000 | 5,622 | 5,834 | 0,964 | 7,475 | 0,752 | 10,764 | 0,522 |
| 10000 | 7,512 | 6,990 | 1,075 | 8,597 | 0,874 | 14,095 | 0,533 |
Рис. 10.7. Графики зависимости ускорения параллельного алгоритма Прима от числа используемых процессоров при различном количестве вершин в модели
Сравнение времени выполнения эксперимента и теоретической оценки Tp из (10.8) приведено в таблице 10.4 и на рис. 10.8.
| 1000 | 0,044 | 0,405 | 0,248 | 0,804 | 0,932 | 1,205 | 1,574 |
| 2000 | 0,208 | 0,820 | 0,684 | 1,613 | 1,800 | 2,412 | 2,159 |
| 3000 | 0,485 | 1,245 | 1,403 | 2,426 | 2,214 | 3,622 | 3,195 |
| 4000 | 0,873 | 1,679 | 1,946 | 3,245 | 3,324 | 4,834 | 5,431 |
| 5000 | 1,432 | 2,122 | 2,665 | 4,068 | 2,933 | 6,048 | 4,119 |
| 6000 | 2,189 | 2,575 | 2,900 | 4,896 | 4,291 | 7,265 | 7,737 |
| 7000 | 3,042 | 3,038 | 3,236 | 5,728 | 6,327 | 8,484 | 8,825 |
| 8000 | 4,150 | 3,510 | 4,462 | 6,566 | 6,993 | 9,705 | 10,390 |
| 9000 | 5,622 | 3,991 | 5,834 | 7,408 | 7,475 | 10,929 | 10,764 |
| 10000 | 7,512 | 4,482 | 6,990 | 8,255 | 8,597 | 12,155 | 14,095 |
Рис. 10.8. Графики экспериментально установленного времени работы параллельного алгоритма Прима и теоретической оценки в зависимости от количества вершин в модели при использовании двух процессоров
Как можно заметить из табл. 10.4 и рис. 10.8, теоретические оценки определяют время выполнения алгоритма Прима с достаточно высокой погрешностью. Причина такого расхождения может состоять в том, что модель Хокни менее точна при оценке времени передачи сообщений с небольшим объемом передаваемых данных. Для уточнения получаемых оценок необходимым является использование других более точных моделей расчета трудоемкости коммуникационных операций – обсуждение этого вопроса проведено в лекции 3.
Сравнение алгоритмов разбиения графов
Рассмотренные алгоритмы разбиения графов различаются точностью получаемых решений, временем выполнения и возможностями для распараллеливания (под точностью понимается величина близости получаемых при помощи алгоритмов решений к оптимальным вариантам разбиения графов). Выбор наиболее подходящего алгоритма в каждом конкретном случае является достаточно сложной и неочевидной задачей. Проведению такого выбора может содействовать сведенная воедино в табл. 10.5 (см. [[67]]) общая характеристика ряда алгоритмов разделения графов, рассмотренных в данном разделе. Дополнительная информация по проблеме оптимального разбиения графов может быть получена, например, в [[67]].
| Покоординатное разбиение | Да | • | • | • • • | |
| Рекурсивный инерционный метод деления пополам | Да | • • | • | • • • | |
| Деление с учетом связности | Нет | • • | • • | • • | |
| Алгоритм Кернигана - Лина | 1 итерация | Нет | • • | • • | • |
| 10 итераций | Нет | • • • | • • • | • • | |
| 50 итераций | Нет | • • • • | • • • • | • • • • | |
Столбец "Необходимость координатной информации" отмечает использование алгоритмом координатной информации об элементах сети или вершинах графа.
Столбец "Точность" дает качественную характеристику величины приближения получаемых алгоритмом решений к оптимальным вариантам разбиения графов. Каждый дополнительный закрашенный кружок определяет примерно 10%-процентное улучшение точности приближения.
Столбец "Время выполнения" показывает относительное время, затрачиваемое различными алгоритмами разбиения. Каждый дополнительный закрашенный кружок соответствует увеличению времени разбиения примерно в 10 раз.
Столбец "Возможности для распараллеливания" характеризует свойства алгоритмов для параллельного выполнения. Алгоритм Кернигана – Лина при выполнении только одной итерации почти не поддается распараллеливанию. Этот же алгоритм при большем количестве итераций, а также метод деления с учетом связности могут быть распараллелены со средней эффективностью. Алгоритм покоординатного разбиения и рекурсивный инерционный метод деления пополам обладают высокими показателями для распараллеливания.
Выделение информационных зависимостей
Выполнение вычислений в подзадачах становится возможным только тогда, когда каждая подзадача (i, j) содержит необходимые для расчетов элементы Aij, Aik, Akj матрицы A. Для исключения дублирования данных разместим в подзадаче (i, j) единственный элемент Aij, тогда получение всех остальных необходимых значений может быть обеспечено только при помощи передачи данных. Таким образом, каждый элемент Akj строки k матрицы A должен быть передан всем подзадачам (k, j), 1jn, а каждый элемент Aik
столбца k матрицы A должен быть передан всем подзадачам (i, k), 1in,– см. рис. 10.3.
Рис. 10.3. Информационная зависимость базовых подзадач (стрелками показаны направления обмена значениями на итерации k)
С учетом выбора базовых подзадач общая схема параллельного выполнения алгоритма Прима будет состоять в следующем:
определяется вершина t графа G, имеющая дугу минимального веса до множества VT. Для выбора такой вершины необходимо осуществить поиск минимума в наборах величин di, имеющихся на каждом из процессоров, и выполнить сборку полученных значений на одном из процессоров;номер выбранной вершины для включения в охватывающее дерево передается всем процессорам;обновляются наборы величин di с учетом добавления новой вершины.
Таким образом, в ходе параллельных вычислений между процессорами выполняются два типа информационных взаимодействий: сбор данных от всех процессоров на одном из процессоров и передача сообщений от одного процессора всем процессорам вычислительной системы.
Задача нахождения минимального охватывающего дерева
Охватывающим деревом (или остовом) неориентированного графа G называется подграф T графа G, который является деревом и содержит все вершины из G. Определив вес подграфа для взвешенного графа как сумму весов входящих в подграф дуг, под минимально охватывающим деревом (МОД) T будем понимать охватывающее дерево минимального веса. Содержательная интерпретация задачи нахождения МОД может состоять, например, в практическом примере построения локальной сети персональных компьютеров с прокладыванием соединительных линий связи минимальной длины. Пример взвешенного неориентированного графа и соответствующего ему минимально охватывающего дерева приведен на рис. 10.6.
Рис. 10.6. Пример взвешенного неориентированного графа (а) и соответствующему ему минимально охватывающего дерева (б)
Дадим общее описание алгоритма решения поставленной задачи, известного под названием метода Прима (the Prim method); более полная информация может быть получена, например, в [[26]].
Задача оптимального разделения графов
Проблема оптимального разделения графов относится к числу часто возникающих задач при проведении различных научных исследований, использующих параллельные вычисления. В качестве примера можно привести задачи обработки данных, в которых области расчетов аппроксимируются двумерными или трехмерными вычислительными сетками. Получение результатов в таких задачах сводится, как правило, к выполнению тех или иных процедур обработки для каждого элемента (узла) сети. При этом в ходе вычислений между соседними элементами сети может происходить передача результатов обработки и т.п. Эффективное решение таких задач на многопроцессорных системах с распределенной памятью предполагает разделение сети между процессорами таким образом, чтобы каждому из процессоров выделялось примерно равное число элементов сети, а межпроцессорные коммуникации, необходимые для выполнения информационного обмена между соседними элементами, были минимальными. На рис. 10.9 показан пример нерегулярной сети, разделенной на 4 части (различные части разбиения сети выделены темным цветом различной интенсивности).
Рис. 10.9. Пример разделения нерегулярной сети
Очевидно, что такие задачи разделения сети между процессорами могут быть сведены к проблеме оптимального разделения графа. Данный подход целесообразен, потому что представление модели вычислений в виде графа позволяет легче решить вопросы хранения обрабатываемых данных и предоставляет возможность применения типовых алгоритмов обработки графов.
Для представления сети в виде графа каждому элементу сети можно поставить в соответствие вершину графа, а дуги графа использовать для отражения свойства близости элементов сети (например, определять дуги между вершинами графа тогда и только тогда, когда соответствующие элементы исходной сети являются соседними). При таком подходе, например, для сети на рис. 10.9, будет сформирован граф, приведенный на рис. 10.10.
Рис. 10.10. Пример графа, моделирующего структуру сети на рис. 10.9
Дополнительная информация по проблеме разделения графов может быть получена, например, в [[67]].
Задача оптимального разделения графов сама может являться предметом распараллеливания. Это бывает необходимо в тех случаях, когда вычислительной мощности и объема оперативной памяти обычных компьютеров недостаточно для эффективного решения задачи. Параллельные алгоритмы разделения графов рассматриваются во многих научных работах: [[20], [38], [44], [48], [49], [65], [74]].
Задача поиска всех кратчайших путей
Исходной информацией для задачи является взвешенный граф G=(V,R), содержащий n вершин (|V|=n), в котором каждому ребру графа приписан неотрицательный вес. Граф будем полагать ориентированным, т.е. если из вершины i есть ребро в вершину j, то из этого не следует наличие ребра из j в i. В случае если вершины все же соединены взаимообратными ребрами, веса, приписанные им, могут не совпадать. Рассмотрим задачу, в которой для имеющегося графа G требуется найти минимальные длины путей между каждой парой вершин графа. В качестве практического примера можно привести задачу составления маршрута движения транспорта между различными городами при заданном расстоянии между населенными пунктами и другие подобные задачи.
В качестве метода, решающего задачу поиска кратчайших путей между всеми парами пунктов назначения, далее используется алгоритм Флойда (the Floyd algorithm) (см, например, [[26]]).
Задачи и упражнения
Используя приведенный программный код, выполните реализацию параллельного алгоритма Флойда. Проведите вычислительные эксперименты. Постройте теоретические оценки с учетом параметров используемой вычислительной системы. Сравните полученные оценки с экспериментальными данными.Выполните реализацию параллельного алгоритма Прима. Проведите вычислительные эксперименты. Постройте теоретические оценки с учетом параметров используемой вычислительной системы. Сравните полученные оценки с экспериментальными данными.Разработайте программную реализацию алгоритма Кернигана – Лина. Дайте оценку возможности распараллеливания этого алгоритма.
Балансировка вычислительной нагрузки процессоров
Как уже отмечалось ранее, вычислительная нагрузка при волновой обработке данных изменяется динамически в ходе вычислений. Данный момент следует учитывать при распределении вычислительной нагрузки между процессорами. Так, например, при фронте волны из 5 блоков и при использовании 4 процессоров обработка волны потребует двух параллельных итераций, во время второй из которых будет задействован только один процессор, а все остальные процессоры будут простаивать, дожидаясь завершения вычислений. Достигнутое ускорение расчетов в этом случае окажется равным 2,5 вместо потенциально возможного значения 4. Подобное снижение эффективности использования процессоров становится менее заметным при большой длине волны, размер которой может регулироваться размером блока. Как общий результат, можно отметить, что размер блока определяет степень разбиения (granularity) вычислений для распараллеливания и является параметром, подбирая значения которого можно управлять эффективностью параллельных вычислений.
Для обеспечения равномерности (балансировки) загрузки процессоров можно задействовать еще один подход, широко используемый для организации параллельных вычислений. Этот подход состоит в том, что все готовые к выполнению в системе вычислительные действия организуются в виде очереди заданий. В ходе вычислений освободившийся процессор может запросить для себя работу из этой очереди; появляющиеся по мере обработки данных новые вычислительные задания пополняют очередь заданий. Такая схема балансировки вычислительной нагрузки между процессорами является простой, наглядной и эффективной. Это позволяет говорить об использовании очереди заданий как общей модели организации параллельных вычислений для систем с общей памятью.
Указанная схема балансировки может быть применена и для рассматриваемого примера. В самом деле, в ходе обработки фронта текущей волны происходит постепенное формирование блоков следующей волны вычислений. Эти блоки могут быть задействованы для обработки при нехватке достаточной вычислительной нагрузки для процессоров.
Общая схема вычислений с использованием очереди заданий может быть представлена в следующем виде:
Алгоритм 11.7. Общая схема вычислений с использованием очереди
// Алгоритм 11.7 // <инициализация служебных данных> // <загрузка в очередь указателя на начальный блок> // взять блок из очереди (если очередь не пуста) while ( (pBlock=GetBlock()) != NULL ) { // <обработка блока> // отметка готовности соседних блоков omp_set_lock(pBlock->pNext.Lock); // сосед справа pBlock->pNext.Count++; if ( pBlock->pNext.Count == 2 ) PutBlock(pBlock->pNext); omp_unset_lock(pBlock->pNext.Lock); omp_set_lock(pBlock->pDown.Lock); // сосед снизу pBlock->pDown.Count++; if ( pBlock->pDown.Count == 2 ) PutBlock(pBlock->pDown); omp_unset_lock(pBlock->pDown.Lock); } // завершение вычислений, т.к. очередь пуста
Для описания имеющихся в задаче блоков узлов сетки в алгоритме используется структура со следующим набором параметров:
Lock – семафор, синхронизирующий доступ к описанию блока;pNext – указатель на соседний справа блок;pDown – указатель на соседний снизу блок;Count – счетчик готовности блока к вычислениям (количество готовых границ блока).
Операции для выборки из очереди и вставки в очередь указателя на готовый к обработке блок узлов сетки обеспечивают соответственно функции GetBlock и PutBlock.
Как следует из приведенной схемы, процессор извлекает блок для обработки из очереди, выполняет необходимые вычисления для блока и отмечает готовность своих границ для соседних справа и снизу блоков. Если при этом оказывается, что у соседних блоков являются подготовленными обе границы, процессор передает эти блоки для запоминания в очередь заданий.
Использование очереди заданий позволяет решить практически все оставшиеся вопросы организации параллельных вычислений для систем с общей памятью. Развитие рассмотренного подхода может предусматривать уточнение правил выделения заданий из очереди для согласования с состояниями процессоров (близкие блоки целесообразно обрабатывать на одних и тех же процессорах), расширение числа имеющихся очередей заданий и т.п.Дополнительная информация по этим вопросам может быть получена, например, в [[63], [76]].
Блочная схема разделения данных
Ленточная схема разделения данных может быть естественным образом обобщена на блочный способ представления сетки области расчетов (см. рис. 11.9). При этом столь радикальное изменение способа разбиения сетки практически не потребует каких-либо существенных корректировок рассмотренной схемы параллельных вычислений. Основной новый момент при блочном представлении данных состоит в увеличении количества граничных строк на каждом процессоре (для блока их количество становится равным 4), что приводит, соответственно, к большему числу операций передачи данных при обмене граничных строк. Сравнивая затраты на организацию передачи граничных строк, можно отметить, что при ленточной схеме для каждого процессора выполняется 4 операции приема-передачи данных, в каждой из которых пересылаются (N+2) значения. Для блочного же способа происходит 8 операций пересылки и объем каждого сообщения равен (N – количество внутренних узлов сетки, NP – число процессоров, размер всех блоков предполагается одинаковым). Тем самым, блочная схема представления области расчетов становится оправданной при большом количестве узлов сетки, когда увеличение количества коммуникационных операций приводит к снижению затрат на пересылку данных в силу сокращения размеров передаваемых сообщений. Результаты экспериментов при блочной схеме разделения данных приведены в табл. 11.5.
| 100 | 210 | 0,06 | 210 | 0,71 | 0,08 | 210 | 0,60 | 0,10 |
| 200 | 273 | 0,35 | 273 | 0,74 | 0,47 | 273 | 1,06 | 0,33 |
| 300 | 305 | 0,92 | 305 | 1,04 | 0,88 | 305 | 2,01 | 0,46 |
| 400 | 318 | 1,69 | 318 | 1,44 | 1,18 | 318 | 2,63 | 0,64 |
| 500 | 343 | 2,88 | 343 | 1,91 | 1,51 | 343 | 3,60 | 0,80 |
| 600 | 336 | 4,04 | 336 | 2,39 | 1,69 | 336 | 4,63 | 0,87 |
| 700 | 344 | 5,68 | 344 | 2,96 | 1,92 | 344 | 5,81 | 0,98 |
| 800 | 343 | 7,37 | 343 | 3,58 | 2,06 | 343 | 7,65 | 0,96 |
| 900 | 358 | 9,94 | 358 | 4,50 | 2,21 | 358 | 9,57 | 1,04 |
| 1000 | 351 | 11,87 | 351 | 4,90 | 2,42 | 351 | 11,16 | 1,06 |
| 2000 | 367 | 50,19 | 367 | 16,07 | 3,12 | 367 | 39,49 | 1,27 |
| 3000 | 364 | 113,17 | 364 | 39,25 | 2,88 | 364 | 85,72 | 1,32 |
(k – количество итераций, t – время (сек), S – ускорение)
При блочном представлении сетки может быть реализован также и волновой метод выполнения расчетов (см. рис. 11.13). Пусть процессоры образуют прямоугольную решетку размером
и процессоры пронумерованы от 0 слева направо по строкам решетки.
Общая схема параллельных вычислений в этом случае имеет вид:
Алгоритм 11.9. Блочная схема разделения данных
Пример 11.9.
(html, txt)
При реализации алгоритма необходимо обеспечить, чтобы в начальный момент времени все процессоры (кроме процессора с нулевым номером) оказались в состоянии передачи своих граничных узлов (верхней строки и левого столбца). Вычисления должен начинать процессор с левым верхним блоком, после завершения обработки которого обновленные значения правого столбца и нижней строки блока нужно переправить правому и нижнему процессорам решетки соответственно. Данные действия обеспечат снятие блокировки процессоров второй диагонали процессорной решетки (ситуация слева на рис. 11.13) и т.д.
Анализ эффективности организации волновых вычислений в системах с распределенной памятью (см. табл. 11.5) показывает значительное снижение полезной вычислительной нагрузки для процессоров, которые занимаются обработкой данных только в моменты, когда их блоки попадают во фронт волны вычислений. При этом балансировка (перераспределение) нагрузки является крайне затруднительной, поскольку связана с пересылкой между процессорами блоков данных большого объема. Возможный интересный способ улучшения ситуации состоит в организации множественной волны вычислений, в соответствии с которой процессоры после отработки волны текущей итерации расчетов могут приступить к выполнению волны следующей итерации метода сеток. Так, например, процессор 0 (см. рис. 11.13), передав после обработки своего блока граничные данные и запустив тем самым вычисления на процессорах 1 и 4, оказывается готовым к исполнению следующей итерации метода Гаусса – Зейделя. После обработки блоков первой (процессоры 1 и 4) и второй (процессор 0) волн к вычислениям окажутся готовыми следующие группы процессоров (для первой волны — процессоры 2, 5 и 8, для второй — процессоры 1 и 4).Кроме того, процессор 0 опять окажется готовым к запуску очередной волны обработки данных. После выполнения NB подобных шагов в обработке будет находиться одновременно NB итераций и все процессоры окажутся задействованными. Подобная схема организации расчетов позволяет рассматривать имеющуюся процессорную решетку как вычислительный конвейер поэтапного выполнения итераций метода сеток. Остановка конвейера может осуществляться, как и ранее, по максимальной погрешности вычислений (проверку условия остановки следует начинать только при достижении полной загрузки конвейера после запуска NB итераций расчетов). Необходимо отметить также, что получаемое после выполнения условия остановки решение задачи Дирихле будет содержать значения узлов сетки от разных итераций метода и не будет, тем самым, совпадать с решением, получаемым при помощи исходного последовательного алгоритма.
Рис. 11.13. Организация волны вычислений при блочной схеме разделения данных
Исключение неоднозначности вычислений
Подход, рассмотренный в п. 11.2.4, уменьшает эффект состязания потоков, но не гарантирует единственности решения при повторении вычислений. Для достижения однозначности необходимо использование дополнительных вычислительных схем.
Возможный и широко применяемый в практике расчетов способ состоит в разделении места хранения результатов вычислений на предыдущей и текущей итерациях метода сеток. Схема такого подхода может быть представлена в следующем общем виде:
Алгоритм 11.4. Параллельная реализация сеточного метода Гаусса – Якоби
// Алгоритм 11.4 omp_lock_t dmax_lock; omp_init_lock(dmax_lock); do { dmax = 0; // максимальное изменение значений u #pragma omp parallel for shared(u,un,N,dmax) private(i,temp,d,dm) for ( i=1; i<N+1; i++ ) { dm = 0; for ( j=1; j<N+1; j++ ) { temp = u[i][j]; un[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-un[i][j]) if ( dm < d ) dm = d; } omp_set_lock(dmax_lock); if ( dmax < dm ) dmax = dm; omp_unset_lock(dmax_lock); } } // конец параллельной области for ( i=1; i<N+1; i++ ) // обновление данных for ( j=1; j<N+1; j++ ) u[i][j] = un[i][j]; } while ( dmax > eps );
Как следует из приведенного алгоритма, результаты предыдущей итерации запоминаются в массиве u, новые вычисления значения запоминаются в дополнительном массиве un. Как результат, независимо от порядка выполнения вычислений для проведения расчетов всегда используются значения величин uij от предыдущей итерации метода. Такая схема реализации сеточных алгоритмов обычно именуется методом Гаусса – Якоби. Этот метод гарантирует однозначность результатов независимо от способа распараллеливания, но требует использования большого дополнительного объема памяти и обладает меньшей (по сравнению с алгоритмом Гаусса – Зейделя) скоростью сходимости. Результаты расчетов с последовательным и параллельным вариантами метода приведены в табл. 11.2.
| 100 | 5257 | 1,39 | 5257 | 0,73 | 1,90 |
| 200 | 23067 | 23,84 | 23067 | 11,00 | 2,17 |
| 300 | 26961 | 226,23 | 26961 | 29,00 | 7,80 |
| 400 | 34377 | 562,94 | 34377 | 66,25 | 8,50 |
| 500 | 56941 | 1330,39 | 56941 | 191,95 | 6,93 |
| 600 | 114342 | 3815,36 | 114342 | 2247,95 | 1,70 |
| 700 | 64433 | 2927,88 | 64433 | 1699,19 | 1,72 |
| 800 | 87099 | 5467,64 | 87099 | 2751,73 | 1,99 |
| 900 | 286188 | 22759,36 | 286188 | 11776,09 | 1,93 |
| 1000 | 152657 | 14258,38 | 152657 | 7397,60 | 1,93 |
| 2000 | 337809 | 134140,64 | 337809 | 70312,45 | 1,91 |
| 3000 | 655210 | 247726,69 | 655210 | 129752,13 | 1,91 |
(k – количество итераций, t – время (сек), S – ускорение)
Рис. . Схема чередования обработки четных и нечетных строк
Иной возможный подход для устранения взаимозависимости параллельных потоков состоит в применении схемы чередования обработки четных и нечетных строк (red/black row alternation scheme), когда выполнение итерации метода сеток подразделяется на два последовательных этапа, на первом из которых обрабатываются строки только с четными номерами, а затем на втором этапе — строки с нечетными номерами (см. рис. 11.7). Данная схема может быть обобщена на применение одновременно и к строкам, и к столбцам (шахматное разбиение) области расчетов.
Рассмотренная схема чередования строк не требует, по сравнению с методом Якоби, какой-либо дополнительной памяти и обеспечивает однозначность решения при многократных запусках программы. Но следует заметить, что оба рассмотренных в данном пункте подхода могут получать результаты, которые не совпадут с решением задачи Дирихле, найденным при помощи последовательного алгоритма. Кроме того, эти вычислительные схемы имеют меньшую область и худшую скорость сходимости, чем исходный вариант метода Гаусса – Зейделя.
Использование OpenMP для организации параллелизма
Рассмотрим возможные способы организации параллельных вычислений для сеточных методов на многопроцессорных вычислительных системах с общей памятью. При изложении материала будем предполагать, что имеющиеся в составе системы процессоры обладают равной производительностью, являются равноправными при доступе к общей памяти и время доступа к памяти является одинаковым (при одновременном доступе нескольких процессоров к одному и тому же элементу памяти очередность и синхронизация доступа обеспечиваются на аппаратном уровне). Как уже отмечалось ранее, многопроцессорные системы подобного типа обычно именуются симметричными мультипроцессорами (symmetric multiprocessors, SMP) – см. п. 1.3.1.
Обычный подход при организации вычислений для подобных систем – создание новых параллельных версий на основе обычных последовательных программ, в которых или автоматически компилятором, или непосредственно программистом выделяются участки независимых друг от друга вычислений. Возможности автоматического анализа программ для порождения параллельных вычислений достаточно ограничены, и второй подход является преобладающим. При этом для разработки параллельных программ могут применяться как новые алгоритмические языки, ориентированные на параллельное программирование, так и уже имеющиеся языки программирования, расширенные некоторым набором операторов для параллельных вычислений.
Оба указанных подхода приводят к необходимости значительной переработки существующего программного обеспечения, и это в значительной степени затрудняет широкое распространение параллельных вычислений. Как результат, в последнее время активно развивается еще один подход к разработке параллельных программ, когда указания программиста по организации параллельных вычислений добавляются в программу при помощи тех или иных внеязыковых средств языка программирования – например, в виде директив или комментариев, которые обрабатываются специальным препроцессором до начала компиляции программы. При этом исходный операторный текст программы остается неизменным, и по нему в случае отсутствия препроцессора компилятор построит исходный последовательный программный код.
Препроцессор же, будучи примененным, заменяет директивы параллелизма на некоторый дополнительный программный код (как правило, в виде обращений к процедурам какой-либо параллельной библиотеки).
Рассмотренный выше подход является основой технологии OpenMP (см., например, [[27]]), наиболее широко применяемой в настоящее время для организации параллельных вычислений на многопроцессорных системах с общей памятью. В рамках данной технологии директивы параллелизма используются для выделения в программе параллельных областей (parallel regions), в которых последовательный исполняемый код может быть разделен на несколько раздельных командных потоков (threads). Далее эти потоки могут исполняться на разных процессорах вычислительной системы. В результате такого подхода программа представляется в виде набора последовательных (однопотоковых) и параллельных (многопотоковых) участков программного кода (см. рис. 11.3). Подобный принцип организации параллелизма получил наименование "вилочного" (fork-join) или пульсирующего параллелизма. Более полная информация по технологии OpenMP может быть получена в литературе (см., например, [[27], [66]]) или в информационных ресурсах сети Интернет. В данной лекции возможности OpenMP будут излагаться в объеме, необходимом для демонстрации возможных способов разработки параллельных программ для рассматриваемого учебного примера решения задачи Дирихле.
Коллективные операции обмена информацией
Для завершения круга вопросов, связанных с параллельной реализацией метода сеток на системах с распределенной памятью, осталось рассмотреть способы вычисления общей для всех процессоров погрешности вычислений. Возможный очевидный подход состоит в передаче всех локальных оценок погрешности, полученных на отдельных полосах сетки, на один какой-либо процессор, вычислении на нем максимального значения и последующей рассылке полученного значения всем процессорам системы. Однако такая схема является крайне неэффективной – количество необходимых операций передачи данных определяется числом процессоров и выполнение этих операций может происходить только в последовательном режиме. Между тем, как показывает анализ требуемых коммуникационных действий, выполнение операций сборки и рассылки данных может быть реализовано с использованием рассмотренной в п. 2.5.2 каскадной схемы обработки данных. На самом деле, получение максимального значения локальных погрешностей, вычисленных на каждом процессоре, может быть обеспечено, например, путем предварительного нахождения максимальных значений для отдельных пар процессоров (такие вычисления могут выполняться параллельно), затем может быть снова осуществлен попарный поиск максимума среди полученных результатов и т.д. Всего, как полагается по каскадной схеме, необходимо выполнить log2NP параллельных итераций для получения конечного значения (NP – количество процессоров).
С учетом возможности применения каскадной схемы для выполнения коллективных операций передачи данных большинство базовых библиотек параллельных программ содержит процедуры для поддержки подобных действий. Так, в стандарте MPI (см. лекцию 5) предусмотрены операции:
Reduce(dm,dmax,op,proc) – процедура сборки на процессоре proc итогового результата dmax среди локальных на каждом процессоре значений dm с применением операции op;Broadcast(dmax,proc) – процедура рассылки с процессора proc значения dmax всем имеющимся процессорам системы.
С учетом перечисленных процедур общая схема вычислений на каждом процессоре может быть представлена в следующем виде:
// Алгоритм 11.8 – уточненный вариант // схема Гаусса-Зейделя, ленточное разделение данных // действия, выполняемые на каждом процессоре do { // обмен граничных строк полос с соседями Sendrecv(u[M][*],N+2,NextProc,u[0][*],N+2,PrevProc); Sendrecv(u[1][*],N+2,PrevProc,u[M+1][*],N+2,NextProc); // <обработка полосы с оценкой погрешности dm> // вычисление общей погрешности вычислений dmax Reduce(dm,dmax,MAX,0); Broadcast(dmax,0); } while ( dmax > eps ); // eps — точность решения
(в приведенном алгоритме переменная dm представляет собой локальную погрешность вычислений на отдельном процессоре, а параметр MAX задает операцию поиска максимального значения для операции сборки). Следует отметить, что в составе MPI имеется процедура Allreduce, которая совмещает действия редукции и рассылки данных. Результаты экспериментов для данного варианта параллельных вычислений для метода Гаусса – Зейделя приведены в рис. 11.4.
Контрольные вопросы
Как определяется задача Дирихле для уравнения Пуассона?В чем состоят основные положения метода конечных разностей?Какие способы распараллеливания сеточных методов могут быть использованы для многопроцессорных вычислительных систем с общей памятью?В каких ситуациях необходима синхронизация параллельных вычислений?Как характеризуется поведение параллельных участков программы при наличии условий состязания потоков?В чем состоит проблема взаимоблокировки?Какие методы могут быть использованы для достижения однозначности результатов параллельных вычислений для сеточных методов?Как изменяется объем вычислений при применении методов волновой обработки данных?Как повысить эффективность методов волновой обработки данных?Как очередь заданий позволяет улучшить балансировку вычислительной нагрузки процессоров?Какие проблемы приходится решать при организации параллельных вычислений на системах с распределенной памятью?Какие основные схемы распределения данных между процессорами могут быть использованы для сеточных методов?Какие основные операции передачи данных используются в параллельных методах решения задачи Дирихле?Каким образом организация множественной волны вычислений позволяет повысить эффективность волновых вычислений в системах с распределенной памятью?
Краткий обзор лекции
В лекции рассматриваются вопросы организации параллельных вычислений для решения задач, в которых при математическом моделировании используются дифференциальные уравнения в частных производных. Для численного решения подобных задач обычно применяется метод конечных разностей (метод сеток), обладающий высокой вычислительной трудоемкостью. В лекции последовательно разбираются возможные способы распараллеливания сеточных методов на многопроцессорных вычислительных системах с общей и распределенной памятью. При этом большое внимание уделяется проблемам, возникающим при организации параллельных вычислений, анализу причин появления таких проблем и нахождению путей их преодоления. Для наглядной демонстрации излагаемого материала в качестве учебного примера рассматривается проблема численного решения задачи Дирихле для уравнения Пуассона.
В пункте 11.1 приводится краткое описание сеточных методов на примере решения задачи Дирихле.
В пункте 11.2 даются возможные способы организации параллельных вычислений при численном решении дифференциальных уравнений в частных производных для вычислительных систем с общей памятью. В основе излагаемого подхода – технология OpenMP, широко применяемая в настоящее время для разработки параллельных программ. В рамках этой технологии параллельный программный код формируется программистом посредством добавления специальных директив или комментариев в существующие последовательные программы. Как результат, программный код является единым для последовательных и параллельных программ, что делает более простым развитие и сопровождение программного обеспечения.
Следует отметить, что принятая в лекции последовательность представления учебного материала может быть рассмотрена как наглядная демонстрация поэтапной методики разработки программного обеспечения. Такой подход позволяет достаточно быстро получать начальные варианты параллельных программ, которые далее могут совершенствоваться для достижения максимально возможной эффективности параллельных вычислений.
В ходе изложения учебного материала в лекции проводится последовательное развитие параллельной программы для решения задачи Дирихле; для каждого очередного варианта программы проводится анализ порождаемых программой параллельных вычислений, определяются причины имеющихся потерь эффективности расчетов и обосновываются пути дальнейшего совершенствования вычислений. Подобный порядок расположения материала позволяет последовательно показать ряд типовых проблем параллельного программирования – излишней синхронизации (serialization), состязания потоков (race condition), тупиков (deadlock) и др. Особое внимание уделяется проблеме возможной неоднозначности результатов последовательных и параллельных вычислений. Для достижения однозначности получаемых результатов расчетов в приводимом учебном материале оценивается возможность применения нескольких различных подходов, последовательный анализ которых приводит к определению методов волновой обработки данных (wavefront or hyperplane methods). На примере реализации волновых схем вычислений дается блочная схема представления данных для эффективного использования быстрой кэш- памяти компьютера. В завершение лекции излагается методика организации очередей заданий для равномерной балансировки вычислительной нагрузки процессоров.
В пункте 11.3 вопросы организации параллельных вычислений при численном решении дифференциальных уравнений в частных производных рассматриваются применительно к вычислительным системам с распределенной памятью. Прежде всего отмечается, что многие проблемы параллельного программирования (состязание вычислений, тупики, сериализация) являются общими для систем с общей и распределенной памятью. Основное отличие параллельных вычислений с распределенной памятью состоит в том, что взаимодействие параллельных участков программы на разных процессорах может быть обеспечено только при помощи передачи сообщений (message passing). При этом эффективность параллельных вычислений во многом определяется равномерностью распределения обрабатываемых данных между процессорами и достигаемой степенью локализации вычислений.
Изложение учебного материала данного раздела лекции начинается с обсуждения общих принципов распределения данных между процессорами, которые применительно к рассматриваемой учебной задаче Дирихле сводятся к одномерной (ленточной) схеме или двумерному (блочному) разбиению области расчетов. Последующее рассмотрение вопросов организации параллельных вычислений проводится в основном на примере ленточной схемы; блочный метод разделения данных представлен в более кратком виде. Среди основных тем, выбранных для обсуждения при изложении ленточной схемы: возможные способы выполнения парных и коллективных операций передачи данных между процессорами, особенности реализации волновых схем вычислений в системах с распределенной памятью, возможность совмещения выполняемых вычислений и операций передачи данных. В завершение лекции проводится сравнительная оценка трудоемкости коммуникационных действий и длительности выполнения вычислительных операций.
Обмен информацией между процессорами
Параллельный вариант метода сеток при ленточном разделении данных состоит в обработке полос на всех имеющихся серверах одновременно в соответствии со следующей схемой работы:
Алгоритм 11.8. Параллельный алгоритм, реализующий метод сеток при ленточном разделении данных
// Алгоритм 11.8 // схема Гаусса-Зейделя, ленточное разделение данных // действия, выполняемые на каждом процессоре do { // <обмен граничных строк полос с соседями> // <обработка полосы> // <вычисление общей погрешности вычислений dmax>} while ( dmax > eps ); // eps — точность решения
Для конкретизации представленных в алгоритме действий введем обозначения:
ProcNum – номер процессора, на котором выполняются описываемые действия,PrevProc, NextProc – номера соседних процессоров, содержащих предшествующую и следующую полосы,NP – количество процессоров,M – количество строк в полосе (без учета продублированных граничных строк),N – количество внутренних узлов в строке сетки (т.е. всего в строке N+2 узла).
При нумерации строк полосы будем считать, что строки 0 и M+1 есть продублированные из соседних полос граничные строки, а строки собственной полосы процессора имеют номера от 1 до M.
Рис. 11.11. Схема передачи граничных строк между соседними процессорами
Процедура обмена граничных строк между соседними процессорами может быть разделена на две последовательные операции, во время первой из которых каждый процессор передает свою нижнюю граничную строку следующему процессору и принимает такую же строку от предыдущего процессора (см. рис. 11.11). Вторая часть передачи строк выполняется в обратном направлении: процессоры передают свои верхние граничные строки своим предыдущим соседям и принимают переданные строки от следующих процессоров.
Выполнение подобных операций передачи данных в общем виде может быть представлено следующим образом (для краткости рассмотрим только первую часть процедуры обмена):
// передача нижней граничной строки следующему // процессору и прием передаваемой строки от // предыдущего процессора if ( ProcNum != NP-1 ) Send(u[M][*],N+2,NextProc); if ( ProcNum != 0 ) Receive(u[0][*],N+2,PrevProc);
(для записи процедур приема-передачи используется близкий к стандарту MPI (см. лекцию 5) формат, где первый и второй параметры представляют пересылаемые данные и их объем, а третий параметр определяет адресата (для операции Send) или источник (для операции Receive) пересылки данных).
Для передачи данных могут быть задействованы два различных механизма. При первом из них выполнение программ, инициировавших операцию передачи, приостанавливается до полного завершения всех действий по пересылке данных (т.е. до момента получения процессором-адресатом всех передаваемых ему данных). Операции приема-передачи, реализуемые подобным образом, обычно называются синхронными или блокирующими. Иной подход – асинхронная или неблокирующая передача — может состоять в том, что операции приема-передачи только инициируют процесс пересылки и на этом завершают свое выполнение. В результате программы, не дожидаясь завершения длительных коммуникационных операций, могут продолжать свои вычислительные действия, проверяя по мере необходимости готовность передаваемых данных. Оба эти варианта операций передачи широко используются при организации параллельных вычислений и имеют свои достоинства и свои недостатки. Синхронные процедуры передачи, как правило, более просты для применения и более надежны; неблокирующие операции могут позволить совместить процессы передачи данных и вычислений, но обычно приводят к повышению сложности программирования. С учетом вышесказанного во всех последующих примерах для организации пересылки данных будут использоваться операции приема-передачи блокирующего типа.
Приведенная выше последовательность блокирующих операций приема-передачи данных (вначале Send, затем Receive) приводит к строго последовательной схеме выполнения процесса пересылок строк, т.к. все процессоры одновременно обращаются к операции Send и переходят в режим ожидания. Первым процессором, который окажется готовым к приему пересылаемых данных, будет сервер с номером NP-1. В результате процессор NP-2 выполнит операцию передачи своей граничной строки и перейдет к приему строки от процессора NP-3 и т.д.
Общее количество повторений таких операций равно NP-1. Аналогично происходит выполнение и второй части процедуры пересылки граничных строк перед началом обработки строк (см. рис. 11.11).
Последовательный характер рассмотренных операций пересылок данных определяется выбранным способом очередности выполнения. Изменим этот порядок очередности при помощи чередования приема и передачи для процессоров с четными и нечетными номерами.
// передача нижней граничной строки следующему // процессору и прием передаваемой строки от // предыдущего процессора if ( ProcNum % 2 == 1 ) { // нечетный процессор if ( ProcNum != NP-1 ) Send(u[M][*],N+2,NextProc); if ( ProcNum != 0 ) Receive(u[0][*],N+2,PrevProc); } else { // процессор с четным номером if ( ProcNum != 0 ) Receive(u[0][*],N+2,PrevProc); if ( ProcNum != NP-1 ) Send(u[M][*],N+2,NextProc); }
Данный прием позволяет выполнить все необходимые операции передачи всего за два последовательных шага. На первом шаге все процессоры с нечетными номерами отправляют данные, а процессоры с четными номерами осуществляют прием этих данных. На втором шаге роли процессоров меняются – четные процессоры выполняют Send, нечетные процессоры исполняют операцию приема Receive.
Рассмотренные последовательности операций приема-передачи для взаимодействия соседних процессоров широко используются в практике параллельных вычислений. Как результат, во многих базовых библиотеках параллельных программ имеются процедуры для поддержки подобных действий. Так, в стандарте MPI (см. лекцию 5) предусмотрена операция Sendrecv, с использованием которой предыдущий фрагмент программного кода может быть записан более кратко:
// передача нижней граничной строки следующему // процессору и прием передаваемой строки от // предыдущего процессора Sendrecv(u[M][*],N+2,NextProc,u[0][*],N+2,PrevProc);
Реализация подобной объединенной функции Sendrecv обычно осуществляется таким образом, чтобы обеспечить и корректную работу на крайних процессорах, когда не нужно выполнять одну из операций передачи или приема, и организацию чередования процедур передачи на процессорах для ухода от тупиковых ситуаций, и возможности параллельного выполнения всех необходимых пересылок данных.
Общие принципы распределения данных
Первая проблема, которую приходится решать при организации параллельных вычислений на системах с распределенной памятью, обычно состоит в выборе способа разделения обрабатываемых данных между вычислительными серверами. Успешность такого разделения определяется достигнутой степенью локализации вычислений на серверах (в силу больших временных задержек при передаче сообщений интенсивность взаимодействия серверов должна быть минимальной).
Рис. 11.10. Ленточное разделение области расчетов между процессорами (кружки представляют граничные узлы сетки)
В рассматриваемом учебном примере по решению задачи Дирихле возможны два различных способа разделения данных – одномерная или ленточная схема (см. рис. 11.10) и двумерное или блочное разбиение (см. рис. 11.9) вычислительной сетки. Дальнейшее изложение учебного материала будет проводиться на примере первого подхода; блочная схема будет рассмотрена позднее в более кратком виде.
При ленточном разбиении область расчетов делится на горизонтальные или вертикальные полосы (не уменьшая общности, далее будем рассматривать только горизонтальные полосы). Число полос определяется количеством процессоров, размер полос обычно является одинаковым, узлы горизонтальных границ (первая и последняя строки) включаются в первую и последнюю полосы соответственно. Полосы для обработки распределяются между процессорами.
Основной момент при организации вычислений с подобным разделением данных состоит в том, что на процессор, выполняющий обработку какой-либо полосы, должны быть продублированы граничные строки предшествующей и следующей полос вычислительной сетки (получаемые в результате расширенные полосы показаны на рис. 11.10 справа пунктирными рамками). Продублированные граничные строки полос используются только при проведении расчетов, пересчет же этих строк происходит в полосах своего исходного месторасположения. Тем самым, дублирование граничных строк должно осуществляться перед началом выполнения каждой очередной итерации метода сеток.
Обзор литературы
Дополнительная информация по численным методам решения дифференциальных уравнений в частных производных может быть получена в [[6], [13]]. Рассмотрение вопросов организации при численном решении дифференциальных уравнений в частных производных проводится в [[5], [51], [60], [63]].
При рассмотрении вопросов организации памяти компьютеров могут оказаться полезными работы [[14], [16]].
Технология MPI для разработки параллельных программ рассмотрена в лекции 5.
Более подробное рассмотрение вопросов, связанных с использованием очередей заданий при организации параллельных вычислений, проводится в [[5], [76]].
Оценка трудоемкости операций передачи данных
Время выполнения коммуникационных операций значительно превышает длительность вычислительных команд. Оценка трудоемкости операций приема-передачи может быть осуществлена с использованием двух основных характеристик сети передачи: латентности (latency), определяющей время подготовки данных к передаче по сети, и пропускной способности сети (bandwidth), задающей объем передаваемых по сети за 1 секунду данных, – более полное изложение вопроса содержится в лекции 3.
Пропускная способность наиболее распространенной на данный момент сети Fast Ethernet – 100 Mбит/с, для более современной сети Gigabit Ethernet – 1000 Мбит/с. В то же время скорость передачи данных в системах с общей памятью обычно составляет сотни и тысячи миллионов байт в секунду. Тем самым, использование систем с распределенной памятью приводит к снижению скорости передачи данных не менее чем в 100 раз.
Еще хуже дело обстоит с латентностью. Для сети Fast Ethernet эта характеристика имеет значение порядка 150 мкс, для сети Gigabit Ethernet – около 100 мкс. Для современных компьютеров с тактовой частотой свыше 2 ГГц различие в производительности достигает не менее чем 10000 – 100000 раз. При указанных характеристиках вычислительной системы для достижения 90% эффективности в рассматриваемом примере решения задачи Дирихле (т.е. чтобы в ходе расчетов обработка данных занимала не менее 90% времени от общей длительности вычислений и только 10% времени тратилось бы на операции передачи данных) размер блоков вычислительной сетки должен быть не менее N=7500 узлов по вертикали и горизонтали (объем вычислений в блоке составляет 5N2
операций с плавающей запятой).
Как результат, можно заключить, что эффективность параллельных вычислений при использовании распределенной памяти определяется в основном интенсивностью и видом выполняемых коммуникационных операций при взаимодействии процессоров. Необходимый при этом анализ параллельных методов и программ может быть выполнен значительно быстрее за счет выделения типовых операций передачи данных – см.
лекцию 3. Так, например, в рассматриваемом учебном примере решения задачи Дирихле практически все пересылки значений сводятся к стандартным коммуникационным действиям, имеющим адекватную поддержку в стандарте MPI (см. рис. 11.14):
Рис. 11.14. Операции передачи данных при выполнении метода сеток с распределенной памятью
рассылка количества узлов сетки всем процессорам – типовая операция передачи данных от одного процессора всем процессорам сети (функция MPI_Bcast);рассылка полос или блоков узлов сетки всем процессорам – типовая операция передачи разных данных от одного процессора всем процессорам сети (функция MPI_Scatter);обмен граничных строк или столбцов сетки между соседними процессорами – типовая операция передачи данных между соседними процессорами сети (функция MPI_Sendrecv);сборка и рассылка погрешности вычислений всем процессорам – типовая операция передачи данных от всех процессоров всем процессорам сети (функция MPI_Allreduce);сборка на одном процессоре решения задачи (всех полос или блоков сетки) – типовая операция передачи данных от всех процессоров сети одному процессору (функция MPI_Gather).
Организация параллельных вычислений для систем с общей памятью
Как следует из приведенного описания, сеточные методы характеризуются значительной вычислительной трудоемкостью
T1=kmN2,
где N есть количество узлов по каждой из координат области D, m — число операций, выполняемых методом для одного узла сетки, k — количество итераций метода до выполнения условия остановки.
Использование процессоров с распределенной памятью является другим общим способом построения многопроцессорных вычислительных систем. Актуальность их становится все более высокой в последнее время в связи с широким развитием высокопроизводительных кластерных вычислительных систем (см. лекцию 1).
Многие проблемы параллельного программирования (состязание вычислений, тупики, сериализация) являются общими для систем с общей и распределенной памятью. Основной момент, который отличает параллельные вычисления с распределенной памятью, состоит в том, что взаимодействие параллельных участков программы на разных процессорах может быть обеспечено только при помощи передачи сообщений (message passing).
Следует отметить, что вычислительный узел системы с распределенной памятью является, как правило, более сложным вычислительным устройством, чем процессор в многопроцессорной системе с общей памятью. Для учета этих различий в дальнейшем процессор с распределенной памятью будет именоваться вычислительным сервером (сервером может быть, в частности, многопроцессорная система с общей памятью). При проведении всех ниже рассмотренных экспериментов использовались 4 компьютера с процессорами Pentium IV, 1300 Mhz, 256 RAM, 100 Mbit Fast Ethernet.
Организация волны вычислений
Представленные в пп. 11.3.1 – 11.3.3 алгоритмы определяют общую схему параллельных вычислений для метода сеток в многопроцессорных системах с распределенной памятью. Далее эта схема может быть конкретизирована реализацией практически всех вариантов методов, рассмотренных для систем с общей памятью (применение дополнительной памяти для схемы Гаусса – Якоби, чередование обработки полос и т.п.). Проработка таких вариантов не привносит каких-либо новых эффектов с точки зрения параллельных вычислений, и их разбор может использоваться как тема заданий для самостоятельных упражнений.
| 100 | 210 | 0,06 | 210 | 0,54 | 0,11 | 210 | 1,27 | 0,05 |
| 200 | 273 | 0,35 | 273 | 0,86 | 0,41 | 273 | 1,37 | 0,26 |
| 300 | 305 | 0,92 | 305 | 0,92 | 1,00 | 305 | 1,83 | 0,50 |
| 400 | 318 | 1,69 | 318 | 1,27 | 1,33 | 318 | 2,53 | 0,67 |
| 500 | 343 | 2,88 | 343 | 1,72 | 1,68 | 343 | 3,26 | 0,88 |
| 600 | 336 | 4,04 | 336 | 2,16 | 1,87 | 336 | 3,66 | 1,10 |
| 700 | 344 | 5,68 | 344 | 2,52 | 2,25 | 344 | 4,64 | 1,22 |
| 800 | 343 | 7,37 | 343 | 3,32 | 2,22 | 343 | 5,65 | 1,30 |
| 900 | 358 | 9,94 | 358 | 4,12 | 2,41 | 358 | 7,53 | 1,32 |
| 1000 | 351 | 11,87 | 351 | 4,43 | 2,68 | 351 | 8,10 | 1,46 |
| 2000 | 367 | 50,19 | 367 | 15,13 | 3,32 | 367 | 27,00 | 1,86 |
| 3000 | 364 | 113,17 | 364 | 37,96 | 2,98 | 364 | 55,76 | 2,03 |
(k – количество итераций, t – время (сек), S – ускорение)
В завершение рассмотрим возможность организации параллельных вычислений, при которых обеспечивалось бы нахождение таких же решений задачи Дирихле, что и при использовании исходного последовательного метода Гаусса – Зейделя. Как отмечалось ранее, такой результат может быть получен за счет организации волновой схемы расчетов. Для образования волны вычислений представим логически каждую полосу узлов области расчетов в виде набора блоков (размер блоков можно положить, в частности, равным ширине полосы) и организуем обработку полос поблочно в последовательном порядке (см. рис. 11.12). Тогда для полного повторения действий последовательного алгоритма вычисления могут быть начаты только для первого блока первой полосы узлов; после того как этот блок будет обработан, для вычислений будут готовы уже два блока – блок 2 первой полосы и блок 1 второй полосы (для обработки блока полосы 2 необходимо передать граничную строку узлов первого блока полосы 1). После обработки указанных блоков к вычислениям будут готовы уже 3 блока, и мы получаем знакомый уже процесс волновой обработки данных (результаты экспериментов см. в табл. 11.4).
Рис. 11.12. Организация волны вычислений при ленточной схеме разделения данных
Интересный момент при организации подобной схемы параллельных вычислений может состоять в попытке совмещения операций пересылки граничных строк и действий по обработке блоков данных.
Параллельные методы решения дифференциальных уравнений в частных производных
Дифференциальные уравнения в частных производных представляют собой широко применяемый математический аппарат при разработке моделей в самых разных областях науки и техники. К сожалению, явное решение этих уравнений в аналитическом виде оказывается возможным только в частных простых случаях, и, как результат, возможность анализа математических моделей, построенных на основе дифференциальных уравнений, обеспечивается при помощи приближенных численных методов решения.
Объем выполняемых при этом вычислений обычно является значительным, и использование высокопроизводительных вычислительных систем традиционно для данной области вычислительной математики.
Проблематика численного решения дифференциальных уравнений в частных производных является областью интенсивных исследований.
Рассмотрим в качестве учебного примера проблему численного решения задачи Дирихле для уравнения Пуассона, которая определяется как задача нахождения функции u=u(x,y), удовлетворяющей в области определения D уравнению
и принимающей значения g(x,y) на границе D0 области D (f и g являются функциями, задаваемыми при постановке задачи). Подобная модель может применяться для описания установившегося течения жидкости, стационарных тепловых полей, процессов теплопередачи с внутренними источниками тепла и деформации упругих пластин. Данный пример часто используется в качестве учебно-практической задачи при изложении возможных способов организации эффективных параллельных вычислений (см. [[12], [60]]).
Для простоты изложения материала в качестве области задания D функции u(x,y) далее будет использоваться единичный квадрат
D={(x,y)D:0x,y1}
Последовательные методы решения задачи Дирихле
Одним из наиболее распространенных подходов к численному решению дифференциальных уравнений является метод конечных разностей (метод сеток) (см., например, [[6], [13], [60]]). Следуя этому подходу, область решения D можно представить в виде дискретного (как правило, равномерного) набора (сетки) точек (узлов). Так, например, прямоугольная сетка в области D может быть задана в виде (рис. 11.1)
где величина N задает количество внутренних узлов по каждой из координат области D.
Обозначим оцениваемую при подобном дискретном представлении аппроксимацию функции u(x,y) в точках (xi, yj) через uij. Тогда, используя пятиточечный шаблон (см. рис. 11.1) для вычисления значений производных, мы можем представить уравнение Пуассона в конечно-разностной форме
Данное уравнение может быть разрешено относительно uij:
uij=0,25(ui-1,j+ui+1,j+ui,j-1-h2fij).
Разностное уравнение, записанное в подобной форме, позволяет определять значение uij по известным значениям функции u(x,y) в соседних узлах используемого шаблона. Данный результат служит основой для построения различных итерационных схем решения задачи Дирихле, в которых в начале вычислений формируется некоторое приближение для значений uij, а затем эти значения последовательно уточняются в соответствии с приведенным соотношением. Так, например, метод Гаусса – Зейделя для проведения итераций уточнения использует правило
по которому очередное k-е приближение значения uij вычисляется по последнему k-му приближению значений ui-1,j и ui,j-1 и предпоследнему (k-1)-му приближению значений ui+1,j и ui,j+1. Выполнение итераций обычно продолжается до тех пор, пока получаемые в результате итераций изменения значений uij не станут меньше некоторой заданной величины (требуемой точности вычислений). Сходимость описанной процедуры (получение решения с любой желаемой точностью) является предметом всестороннего математического анализа (см., например, [[6], [13], [60]]), здесь же отметим, что последовательность решений, получаемых методом сеток, равномерно сходится к решению задачи Дирихле, а погрешность решения имеет порядок h2.
Рис. 11.1. Прямоугольная сетка в области D ( темные точки представляют внутренние узлы сетки, нумерация узлов в строках слева направо, а в столбцах — сверху вниз)
Рассмотренный алгоритм (метод Гаусса – Зейделя) на псевдокоде, приближенном к алгоритмическому языку С++, может быть представлен в виде:
Алгоритм 11.1. Последовательный алгоритм Гаусса – Зейделя
// Алгоритм 11.1 do { dmax = 0; // максимальное изменение значений u for ( i=1; i<N+1; i++ ) for ( j=1; j<N+1; j++ ) { temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); dm = fabs(temp-u[i][j]); if ( dmax < dm ) dmax = dm; } } while ( dmax > eps );
(напомним, что значения uij при индексах i,j=0,N+1 являются граничными, задаются при постановке задачи и не изменяются в ходе вычислений).
Рис. 11.2. Вид функции u(x,y) в примере для задачи Дирихле
Для примера на рис. 11.2
приведен вид функции u(x,y), полученной для задачи Дирихле при следующих граничных условиях:
Общее количество итераций метода Гаусса – Зейделя составило 210 при точности решения eps=0,1 и N=100 (в качестве начального приближения величин uij использовались значения, сгенерированные датчиком случайных чисел из диапазона [-100, 100]).
pragma omp parallel for
| // Алгоритм 11.5 omp_lock_t dmax_lock; omp_init_lock(dmax_lock); do { dmax = 0; // максимальное изменение значений u // нарастание волны (nx – размер волны) for ( nx=1; nx<N+1; nx++ ) { dm[nx] = 0; # pragma omp parallel for shared(u,N,nx,dm) private(i,j,temp,d) for ( i=1; i<nx+1; i++ ) { j = nx + 1 – i; temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-u[i][j]) if ( dm[i] < d ) dm[i] = d; } // конец параллельной области } // затухание волны for ( nx=N-1; nx>0; nx-- ) { #pragma omp parallel for shared(u,N,nx,dm) private(i,j,temp,d) for ( i=N-nx+1; i<N+1; i++ ) { j = 2*N - nx – I + 1; temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-u[i][j]) if ( dm[i] < d ) dm[i] = d; } // конец параллельной области } #pragma omp parallel for shared(nx,dm,dmax) private(i) for ( i=1; i<nx+1; i++ ) { omp_set_lock(dmax_lock); if ( dmax < dm[i] ) dmax = dm[i]; omp_unset_lock(dmax_lock); } // конец параллельной области } while ( dmax > eps ); |
| Пример 11.5. |
| Закрыть окно |
// Алгоритм 11.5
omp_lock_t dmax_lock;
omp_init_lock(dmax_lock);
do {
dmax = 0; // максимальное изменение значений u
// нарастание волны (nx – размер волны)
for ( nx=1; nx
dm[nx] = 0;
# pragma omp parallel for shared(u,N,nx,dm) private(i,j,temp,d)
for ( i=1; i
j = nx + 1 – i;
temp = u[i][j];
u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+
u[i][j-1]+u[i][j+1]–h*h*f[i][j]);
d = fabs(temp-u[i][j])
if ( dm[i] < d ) dm[i] = d;
} // конец параллельной области
}
// затухание волны
for ( nx=N-1; nx>0; nx-- ) {
#pragma omp parallel for shared(u,N,nx,dm) private(i,j,temp,d)
for ( i=N-nx+1; i
j = 2*N - nx – I + 1;
temp = u[i][j];
u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+
u[i][j-1]+u[i][j+1]–h*h*f[i][j]);
d = fabs(temp-u[i][j])
if ( dm[i] < d ) dm[i] = d;
} // конец параллельной области
}
#pragma omp parallel for shared(nx,dm,dmax) private(i)
for ( i=1; i
omp_set_lock(dmax_lock);
if ( dmax < dm[i] ) dmax = dm[i];
omp_unset_lock(dmax_lock);
} // конец параллельной области
} while ( dmax > eps );
Проблема синхронизации параллельных вычислений
Первый вариант параллельного алгоритма для метода сеток может быть получен, если разрешить произвольный порядок пересчета значений uij. Программа для данного способа вычислений может быть представлена в следующем виде:
Алгоритм 11.2. Первый вариант параллельного алгоритма Гаусса – Зейделя
// Алгоритм 11.2 omp_lock_t dmax_lock; omp_init_lock (dmax_lock); do { dmax = 0; // максимальное изменение значений u #pragma omp parallel for shared(u,N,dmax) private(i,temp,d) for ( i=1; i<N+1; i++ ) { #pragma omp parallel for shared(u,N,dmax) private(j,temp,d) for ( j=1; j<N+1; j++ ) { temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-u[i][j]) omp_set_lock(dmax_lock); if ( dmax < d ) dmax = d; omp_unset_lock(dmax_lock); } // конец вложенной параллельной области } // конец внешней параллельной области } while ( dmax > eps );
Следует отметить, что программа получена из исходного последовательного кода путем добавления директив и операторов обращения к функциям библиотеки OpenMP (эти дополнительные строки в программе выделены темным шрифтом).
Рис. 11.3. Параллельные области, создаваемые директивами OpenMP
Как следует из текста программы, параллельные области в данном примере задаются директивой parallel for, являются вложенными и включают в свой состав операторы цикла for. Компилятор, поддерживающий технологию OpenMP, разделяет выполнение итераций цикла между несколькими потоками программы, количество которых обычно совпадает с числом процессоров в вычислительной системе. Параметры директивы shared и private определяют доступность данных в потоках программы – переменные, описанные как shared, являются общими для потоков, для переменных с описанием private создаются отдельные копии для каждого потока, которые могут использоваться в потоках независимо друг от друга.
Наличие общих данных обеспечивает возможность взаимодействия потоков. В этом плане разделяемые переменные могут рассматриваться как общие ресурсы потоков, и, как результат, их применение должно выполняться с соблюдением правил взаимоисключения (изменение каким-либо потоком значений общих переменных должно приводить к блокировке доступа к модифицируемым данным для всех остальных потоков).
В данном примере таким разделяемым ресурсом является величина dmax, доступ потоков к которой регулируется специальной служебной переменной (замком) dmax_lock и функциями omp_set_lock (разрешение или блокировка доступа) и omp_unset_lock (снятие запрета на доступ). Подобная организация программы гарантирует единственность доступа потоков для изменения разделяемых данных. Участки программного кода (блоки между обращениями к функциям omp_set_lock и omp_unset_lock), для которых обеспечивается взаимоисключение, обычно именуются критическими секциями.
Результаты вычислительных экспериментов приведены в табл. 11.1
(здесь и далее для параллельных программ, разработанных с использованием технологии OpenMP, использовался четырехпроцессорный сервер кластера Нижегородского университета с процессорами Pentium III, 700 Mhz, 512 RAM).
Оценим полученный результат. Разработанный параллельный алгоритм является корректным, т.е. обеспечивающим решение поставленной задачи. Использованный при разработке подход обеспечивает достижение практически максимально возможного параллелизма – для выполнения программы может быть задействовано вплоть до N2 процессоров. Тем не менее результат не может быть признан удовлетворительным: программа будет работать медленно и вместо ускорения мы получим замедление вычислений. Основная причина такого положения дел – чрезмерно высокая синхронизация параллельных участков программы. В нашем примере каждый параллельный поток после усреднения значений uij должен проверить (и возможно, изменить) значение величины dmax. Разрешение на использование переменной может получить только один поток – все остальные потоки должны быть блокированы. После освобождения общей переменной управление может получить следующий поток и т.д. В результате необходимости синхронизации доступа многопотоковая параллельная программа превращается фактически в последовательно выполняемый код, причем менее эффективный, чем исходный последовательный вариант, т.к. организация синхронизации приводит к дополнительным вычислительным затратам – см.
рис. 11.4. Следует обратить внимание, что, несмотря на идеальное распределение вычислительной нагрузки между процессорами, для приведенного на рис. 11.4 соотношения параллельных и последовательных вычислений, в каждый текущий момент времени (после момента первой синхронизации) только не более двух процессоров одновременно выполняют действия, связанные с решением задачи. Подобный эффект вырождения параллелизма из-за интенсивной синхронизации параллельных участков программы обычно именуется сериализацией (serialization).
| 100 | 210 | 0,06 | 210 | 1,97 | 0,03 | 210 | 0,03 | 2,03 |
| 200 | 273 | 0,34 | 273 | 11,22 | 0,03 | 273 | 0,14 | 2,43 |
| 300 | 305 | 0,88 | 305 | 29,09 | 0,03 | 305 | 0,36 | 2,43 |
| 400 | 318 | 3,78 | 318 | 54,20 | 0,07 | 318 | 0,64 | 5,90 |
| 500 | 343 | 6,00 | 343 | 85,84 | 0,07 | 343 | 1,06 | 5,64 |
| 600 | 336 | 8,81 | 336 | 126,38 | 0,07 | 336 | 1,50 | 5,88 |
| 700 | 344 | 12,11 | 344 | 178,30 | 0,07 | 344 | 2,42 | 5,00 |
| 800 | 343 | 16,41 | 343 | 234,70 | 0,07 | 343 | 8,08 | 2,03 |
| 900 | 358 | 20,61 | 358 | 295,03 | 0,07 | 358 | 11,03 | 1,87 |
| 1000 | 351 | 25,59 | 351 | 366,16 | 0,07 | 351 | 13,69 | 1,87 |
| 2000 | 367 | 106,75 | 367 | 1585,84 | 0,07 | 367 | 56,63 | 1,89 |
| 3000 | 370 | 243,00 | 370 | 3598,53 | 0,07 | 370 | 128,66 | 1,89 |
(k-количество итераций, t-время (сек), S-ускорение)
Рис. 11.4. Пример возможной схемы выполнения параллельных потоков при наличии синхронизации (взаимоисключения)
Как показывают выполненные рассуждения, путь для достижения эффективности параллельных вычислений лежит в уменьшении необходимых моментов синхронизации параллельных участков программы. Так, в нашем примере мы можем ограничиться распараллеливанием только одного внешнего цикла for. Кроме того, для снижения количества возможных блокировок применим для оценки максимальной погрешности многоуровневую схему расчета: пусть параллельно выполняемый поток первоначально формирует локальную оценку погрешности dm только для своих обрабатываемых данных (одной или нескольких строк сетки), затем при завершении вычислений поток сравнивает свою оценку dm с общей оценкой погрешности dmax.
Новый вариант программы решения задачи Дирихле имеет вид:
Алгоритм 11.3. Второй вариант параллельного алгоритма Гаусса – Зейделя
// Алгоритм 11.3 omp_lock_t dmax_lock; omp_init_lock(dmax_lock); do { dmax = 0; // максимальное изменение значений u #pragma omp parallel for shared(u,N,dmax) private(i,temp,d,dm) for ( i=1; i<N+1; i++ ) { dm = 0; for ( j=1; j<N+1; j++ ) { temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); d = fabs(temp-u[i][j]) if ( dm < d ) dm = d; } omp_set_lock(dmax_lock); if ( dmax < dm ) dmax = dm; omp_unset_lock(dmax_lock); } } // конец параллельной области } while ( dmax > eps );
Как результат выполненного изменения схемы вычислений, количество обращений к общей переменной dmax уменьшается с N2 до N раз, что должно приводить к существенному снижению затрат на синхронизацию потоков и уменьшению проявления эффекта сериализации вычислений. Результаты экспериментов с данным вариантом параллельного алгоритма, приведенные в табл. 11.1, показывают существенное изменение ситуации – ускорение в ряде экспериментов оказывается даже большим, чем используемое количество процессоров (такой эффект сверхлинейного ускорения достигается за счет наличия у каждого из процессоров вычислительного сервера своей быстрой кэш-памяти). Следует также обратить внимание, что улучшение показателей параллельного алгоритма достигнуто при снижении максимально возможного параллелизма (для выполнения программы может использоваться не более N процессоров).
Проблема взаимоблокировки
Возможный подход для получения однозначных результатов (уход от состязания потоков) может состоять в ограничении доступа к узлам сетки, которые обрабатываются в параллельных потоках программы. Для этого можно ввести набор замков row_lock[N], который позволит потокам закрывать доступ к "своим" строкам сетки.
// поток обрабатывает i строку сетки omp_set_lock(row_lock[i]); omp_set_lock(row_lock[i+1]); omp_set_lock(row_lock[i-1]); // обработка i строки сетки omp_unset_lock(row_lock[i]); omp_unset_lock(row_lock[i+1]); omp_unset_lock(row_lock[i-1]);
Закрыв доступ к своим данным, параллельный поток уже не будет зависеть от динамики выполнения других параллельных участков программы. Результат вычислений потока однозначно определяется значениями данных в момент начала расчетов.
Данный подход позволяет продемонстрировать еще одну проблему, которая может возникать в ходе параллельных вычислений. Эта проблема состоит в том, что при организации доступа к множественным общим переменным может появляться конфликт между параллельными потоками и этот конфликт не может быть разрешен успешно. Так, в приведенном фрагменте программного кода при обработке потоками двух последовательных строк (например, строк 1 и 2) может сложиться ситуация, когда потоки блокируют сначала строки 1 и 2 и только затем переходят к блокировке оставшихся строк (см. рис. 11.6). В этом случае доступ к необходимым строкам не может быть обеспечен ни для одного потока – возникает неразрешимая ситуация, обычно именуемая тупиком. Как можно показать, необходимым условием тупика является наличие цикла в графе распределения и запросов ресурсов. В рассматриваемом примере уход от цикла может состоять в строго последовательной схеме блокировки строк потока.
// поток обрабатывает i строку сетки omp_set_lock(row_lock[i+1]); omp_set_lock(row_lock[i]); omp_set_lock(row_lock[i-1]); // <обработка i строки сетки> omp_unset_lock(row_lock[i+1]); omp_unset_lock(row_lock[i]); omp_unset_lock(row_lock[i-1]);
(следует отметить, что и эта схема блокировки строк может оказаться тупиковой, если рассматривать модифицированную задачу Дирихле, в которой горизонтальные границы являются "склеенными").
Рис. 11.6. Ситуация тупика при доступе к строкам сетки (поток 1 владеет строкой 1 и запрашивает строку 2, поток 2 владеет строкой 2 и запрашивает строку 1)
Волновые схемы параллельных вычислений
Рассмотрим теперь возможность построения параллельного алгоритма, который выполнял бы только те вычислительные действия, что и последовательный метод (может быть только в несколько ином порядке) и, как результат, обеспечивал бы получение точно таких же решений исходной вычислительной задачи. Как уже было отмечено выше, в последовательном алгоритме каждое очередное k-е приближение значения ui,j вычисляется по последнему k-му приближению значений ui-1,j и ui,j-1 и предпоследнему (k-1)-му приближению значений ui+1,j и ui,j+1. Таким образом, при требовании совпадения результатов вычислений последовательных и параллельных вычислительных схем в начале каждой итерации метода только одно значение u11 может быть пересчитано (возможности для распараллеливания нет). Но далее после пересчета u11 вычисления могут выполняться уже в двух узлах сетки u12 и u21 (в этих узлах выполняются условия последовательной схемы), затем после пересчета узлов u12 и u21 — в узлах u13, u22 и u31 и т.д. Обобщая сказанное, можно увидеть, что выполнение итерации метода сеток можно разбить на последовательность шагов, на каждом из которых к вычислениям окажутся подготовленными узлы вспомогательной диагонали сетки с номером, определяемым номером этапа – см. рис. 11.8. Получаемая в результате вычислительная схема получила наименование волны или фронта вычислений, а алгоритмы, построенные на ее основе, — методов волновой обработки данных (wavefront или hyperplane methods). Следует отметить, что в нашем случае размер волны (степень возможного параллелизма) динамически изменяется в ходе вычислений – волна нарастает до своего пика, а затем затухает при приближении к правому нижнему узлу сетки.
Рис. 11.8. Движение фронта волны вычислений
Возможная схема параллельного метода, основанного на эффекте волны вычислений, может быть представлена в следующей форме.
Алгоритм 11.5. Параллельный алгоритм, реализующий волновую схему вычислений
Пример 11.5.
(html, txt)
При разработке алгоритма, реализующего волновую схему вычислений, оценку погрешности решения можно осуществлять для каждой строки в отдельности (массив значений dm).
Этот массив является общим для всех выполняемых потоков, однако синхронизации доступа к элементам не требуется, так как потоки используют всегда разные элементы массива (фронт волны вычислений содержит только по одному узлу строк сетки).
После обработки всех элементов волны в составе массива dm находится максимальная погрешность выполненной итерации вычислений. Однако именно эта последняя часть расчетов может оказаться наиболее неэффективной из-за высоких дополнительных затрат на синхронизацию. Улучшение ситуации, как и ранее, может быть достигнуто за счет увеличения размера последовательных участков и сокращения, тем самым, количества необходимых взаимодействий параллельных участков вычислений. Возможный вариант реализации такого подхода может состоять в следующем:
chunk = 200; // размер последовательного участка #pragma omp parallel for shared(n,dm,dmax) private(i,d) for ( i=1; i<nx+1; i+=chunk ) { d = 0; for ( j=i; j<i+chunk; j++ ) if ( d < dm[j] ) d = dm[j]; omp_set_lock(dmax_lock); if ( dmax < d ) dmax = d; omp_unset_lock(dmax_lock); } // конец параллельной области
| 100 | 210 | 0,06 | 210 | 0,30 | 0,21 | 210 | 0,16 | 0,40 |
| 200 | 273 | 0,34 | 273 | 0,86 | 0,40 | 273 | 0,59 | 0,58 |
| 300 | 0,88 | 305 | 1,63 | 0,54 | 305 | 1,53 | 0,57 | |
| 400 | 318 | 3,78 | 318 | 2,50 | 1,51 | 318 | 2,36 | 1,60 |
| 500 | 343 | 6,00 | 343 | 3,53 | 1,70 | 343 | 4,03 | 1,49 |
| 600 | 336 | 8,81 | 336 | 5,20 | 1,69 | 336 | 5,34 | 1,65 |
| 700 | 344 | 12,11 | 344 | 8,13 | 1,49 | 344 | 10,00 | 1,21 |
| 800 | 343 | 16,41 | 343 | 12,08 | 1,36 | 343 | 12,64 | 1,30 |
| 900 | 358 | 20,61 | 358 | 14,98 | 1,38 | 358 | 15,59 | 1,32 |
| 1000 | 351 | 25,59 | 351 | 18,27 | 1,40 | 351 | 19,30 | 1,33 |
| 2000 | 367 | 106,75 | 367 | 69,08 | 1,55 | 367 | 65,72 | 1,62 |
| 3000 | 370 | 243,00 | 370 | 149,36 | 1,63 | 370 | 140,89 | 1,72 |
(k – количество итераций, t – время (сек), S – ускорение)
Подобный прием укрупнения последовательных участков вычислений для снижения затрат на синхронизацию именуется фрагментированием (chunking).
Возможность неоднозначности вычислений в параллельных программах
Последний рассмотренный вариант организации параллельных вычислений для метода сеток обеспечивает практически максимально возможное ускорение выполняемых расчетов – так, в экспериментах данное ускорение достигало величины 5,9 при использовании четырехпроцессорного вычислительного сервера. Вместе с этим необходимо отметить, что разработанная вычислительная схема расчетов имеет важную принципиальную особенность: порождаемая при вычислениях последовательность обработки данных может различаться при разных запусках программы даже при одних и тех же исходных параметрах решаемой задачи. Данный эффект может проявляться в силу изменения каких-либо условий выполнения программы (вычислительной нагрузки, алгоритмов синхронизации потоков и т.п.), что может повлиять на временные соотношения между потоками (см. рис. 11.5). Взаиморасположение потоков по области расчетов может быть различным: одни потоки могут опережать другие и, наоборот, часть потоков могут отставать (при этом характер взаиморасположения может меняться в ходе вычислений). Подобное поведение параллельных участков программы обычно именуется состязанием потоков (race condition) и отражает важный принцип параллельного программирования – временная динамика выполнения параллельных потоков не должна учитываться при разработке параллельных алгоритмов и программ.
Рис. 11.5. Возможные различные варианты взаиморасположения параллельных потоков (состязание потоков)
В рассматриваемом примере при вычислении нового значения uij в зависимости от условий выполнения могут использоваться разные (от предыдущей или текущей итераций) оценки соседних значений по вертикали. Тем самым, количество итераций метода до выполнения условия остановки и, самое главное, конечное решение задачи могут различаться при повторных запусках программы. Получаемые оценки величин uij. будут соответствовать точному решению задачи в пределах задаваемой точности, но, тем не менее, могут быть различными. Применение вычислений такого типа для сеточных алгоритмов получило наименование метода хаотической релаксации (chaotic relaxation).
Задачи и упражнения
Выполните реализацию первого и второго вариантов параллельного алгоритма Гаусса – Зейделя для систем с общей памятью. Проведите вычислительные эксперименты и сравните время выполнения разработанных программ.Выполните реализации параллельного алгоритма Гаусса – Зейделя при волновой схеме организации вычислений и блочном представлении обрабатываемых данных. Проведите вычислительные эксперименты при разном размере блоков и сравните получаемые характеристики эффективности параллельных вычислений.Выполните реализацию очереди заданий для параллельного алгоритма Гаусса – Зейделя. Подготовьте несколько разных правил выделения заданий из очереди и проведите оценку эффективности для каждого использованного правила.
