Участник:Ilya4870/Алгоритм разбиения графа методом рекурсивной координатной бисекции

1 Свойства и структура алгоритмов

1.1 Общее описание алгоритма

Задача рациональной декомпозиции расчетных сеток возникает при численном моделировании на высокопроизводительных вычислительных системах проблем механики сплошных сред, импульсной энергетики, электродинамики и многих других. Такие задачи сводятся к более общей задачи - разбиение вершин графа. Оптимальным считается разбиение на домены, при котором выровнен суммарный вес вершин в доменах и минимизирован суммарный вес разрезанных ребер (разрезанное ребро – ребро, соединяющее вершины из разных доменов). Одним из методов разбиения вершин графа на равные по мощности блоки является метод рекурсивной координатной бисекции.

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

1.2 Математическое описание алгоритма

Пусть [math]G = (V,E)[/math] – неориентированный граф, где [math]V = \{v_i\}[/math] – множество вершин, [math]E = \{e_{ij}\}[/math] – множество рёбер. И вершины, и рёбра, имеют целочисленные веса [math]|v_i|[/math] и [math]|e_{ij}|[/math], соответственно.

Декомпозиция является отображением множества вершин [math]V[/math] на заданное число [math]p[/math] доменов [math]S_j[/math] такое, что каждая вершина [math]v_i[/math] принадлежит некоторому домену [math]S_j[/math], и [math]S_k \cap S_m = \emptyset[/math].

Суммарный вес вершин в домене [math]S_j[/math] равен сумме весов вершин, принадлежащих этому домену: [math]|S_j| = \sum_{i}{|v_i|}, v_i \in S_j[/math].

Суммарный вес разрезанных рёбер равен сумме весов рёбер, соединяющих вершины из разных доменов: [math]|E_c|=\sum_{ij}{|e_ij|}, v_i \in S_k, v_j \in S_m, k \lt m[/math].

Каждый домен [math]S_j[/math] является некоторым подграфом [math]G^j = (V^j, E^j)[/math] графа [math]G = (V, E)[/math], где [math]V^j \subseteq V, E^j \subseteq E[/math]. Домен является связным, если его граф состоит из одной компоненты связности, то есть между любой парой вершин этого графа существует как минимум один путь.

Требуется найти такое разбиение множества вершин [math]V[/math] на заданное число [math]p[/math] связных доменов [math]S_j[/math], при котором:

• выровнен суммарный вес вершин в доменах: [math]|S_j| \approx |V|/p[/math],[math](1)[/math]

• минимизирован суммарный вес разрезанных рёбер: [math]min(|E_c|)[/math].[math](2)[/math]

Используя обозначения выше, определим:

• Входные данные алгоритма: множество координат вершин [math]C = \{c_i\}[/math] графа [math]G = (V, E)[/math], где [math]c_i = (v_i^x, v_i^y, v_i^z)[/math] - координаты вершины [math]v_i \in V[/math], число доменов [math]p[/math].

• Выходные данные алгоритма: разбиение множества вершин [math]V[/math] на заданное число [math]p[/math] доменов [math]S_j[/math].

Основные этапы алгоритма:

• вычисляются протяжённости осей: [math]l_k = |max(v_i^k) - min(v_i^k)|[/math] для [math]k = x, y, z[/math],

• производится сортировка множества [math]C[/math] по координате, соответствующей наибольшей [math]l_k[/math],

• множество [math]C[/math] разбивается на две части в соответствии с новым порядком элементов, то есть перпендикулярно координатной оси, вдоль которой она имеет наибольшую протяжённость; соотношение размеров частей зависит от количества доменов, которые должны быть образованы в каждой из частей,

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

1.3 Вычислительное ядро алгоритма

Большую часть времени занимает сортировка массива координат вершин [math]C[/math], выполняемая на каждом шаге рекурсии. Сортировка может быть выполнена любым способом, в зависимости от реализации. В дальнейшем положим, что используется алгоритм с оптимальной вычислительной сложностью [math]O(nlog_2{n})[/math].

1.4 Макроструктура алгоритма

Как записано и в описании ядра алгоритма, основную часть метода составляют сортировки массивов координат вершин [math]C[/math], выполняемые на каждом шаге рекурсии.

1.5 Схема реализации последовательного алгоритма

1. Вычисляются минимумы и максимумы по каждой из координат, определяющие параллелепипед, охватывающий вершины данной области. Далее запускается процесс рекурсии, на каждом шаге которого выполняется следующее:
   1. Определяется координата [math]j[/math], вдоль которой параллелепипед, обрабатывающийся на данном этапе, имеет наибольшую протяженность.
   2. Определяется количество доменов, которое делится на данном этапе. В левую часть (вершины с меньшей координатой [math]j[/math]) попадет половина доменов. Сначала набираются домены нулевого типа, потом, если они закончились, первого, и т.д. до получения нужного числа доменов в левой части. Подсчитывается общее число вершин n1, которое окажется в этих доменах. В правой части будут сформированы оставшиеся домены.
   3. Пусть все вершины в параллелепипеде имеют координаты [math]j[/math], принадлежащие интервалу [math][j_{min}, j_{max}][/math]. Интервал делится на [math]K[/math] равных малых интервалов. [math]K[/math] – некоторая константа. Подсчитываются количества вершин, принадлежащие каждому малому интервалу. Исходя из того, что в левую часть должна попасть [math]n_1[/math] вершина с меньшей координатой [math]j[/math], определяется, в какой малый интервал [math][j_{from}, j_{to}][/math] попадает медиана. Затем вершины из интервала [math][j_{from}, j_{to}][/math] записываются в отдельный массив, отсортированный по всем координатам. Определяется медиана – последняя вершина, которая попадет в левую часть.
   4. Если предполагается построение локальных деревьев, для вершины [math]i[/math] двоичного дерева запоминаются данные о медиане: ее координаты, номер вершины медианы и номер координаты [math]j[/math], по которой проводилось деление.
   5. Проводится разделение данных по аналогии с обменной сортировкой с разделением Хоару (быстрая сортировка), чтобы на следующих этапах просматривались только вершины из разбиваемой области. Два указателя, [math]begin[/math] и [math]end[/math], выставляются на начало и на конец массива данных, соответственно. Вершина, на которую указывает указатель [math]begin[/math], сравнивается с медианой по всем координатам, и, если она попадает в левую часть параллелепипеда, указатель [math]begin[/math] сдвигается к следующей вершине. Так до тех пор, пока указатель [math]begin[/math] не остановится на вершине, которая должна попасть в правую часть. Указатель [math]end[/math] конца массива сдвигается влево до тех пор, пока не остановится на вершине, которая должна попасть в левую часть. Если указатель [math]begin[/math] находится левее указателя [math]end[/math], вершины, на которые указывают указатели, меняются местами, после чего указатель [math]begin[/math] сдвигается вправо, а [math]end[/math] – влево.Если указатель [math]begin[/math] все еще левее указателя [math]end[/math], весь процесс повторяется снова. Если после выхода из цикла указатели указывают на одну и ту же вершину, она сравнивается с медианой по всем координатам, и если вершина попадает в левую часть, указатель [math]begin[/math] сдвигается к следующей вершине. Указатель [math]begin[/math] определяет начало массива данных правой части.
2. По рекурсии левая и правая части отправляются на деление. Левой части передается указатель на начало исходного массива данных и число вершин [math]n_1[/math]. Максимальные координаты параллелепипеда остаются теми же, только максимальная координата [math]j[/math] приравнивается к координате [math]j[/math] медианы. Передается число доменов каждого вида, которое будет в левой части. Номер вершины двоичного дерева, отвечающей за деление левой части, будет равным [math]2i[/math]. Правой части передается указатель [math]begin[/math] в качестве нового указателя на массив данных и оставшееся число вершин. Минимальная координата [math]j[/math] приравнивается к координате медианы. Передаются числа доменов каждого вида. Номер вершины двоичного дерева, отвечающей за деление правой части, будет равным [math]2i + 1[/math].
3. Перед делением каждой из частей проверяется, сколько доменов осталось в данной части. Если остался только один домен, то все вершины помещаются в данный домен. Домены считаются по порядку от номера [math]gr_0[/math] первого домена. Если предполагается построение локальных деревьев, в лист дерева записывается номер домена. Запоминается максимальный номер вершины. Если в части осталось больше одного домена, она отправляется на деление.

Рисунок 1. Пример координатной бисекции на три домена.

Алгоритм работает только с координатами вершин и не учитывает связи между ними, что делает его экономичным по памяти.

1.6 Последовательная сложность алгоритма

Для разбиения графа с количеством вершин [math]n[/math] методом рекурсивной координатной бисекции на [math]p[/math] доменов требуется:

• для каждого шага рекурсии:
  • определить максимум и минимум координат, [math]O(n)[/math] операций сравнения,
  • отсортировать подобласть по одной из координат, [math]O(nlog_2{n})[/math] операций сравнения.

Глубина рекурсии - [math]log_2{p}+1[/math].

Сравнения составляют основную часть алгоритма.

Таким образом, алгоритм имеет последовательную сложность [math]O(n, p) = n \log_2{n} \log_2{p}.[/math]

1.7 Информационный граф

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

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

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

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

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

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

Рис.1. Информационная структура алгоритма умножения матриц

Рис.2. Информационная структура одного из вариантов алгоритма решения систем линейных алгебраических уравнений с блочно-двухдиагональной матрицей

1.8 Ресурс параллелизма алгоритма

Исходя из того, что метод является геометрическим, алгоритм обладает координатным параллелизмом. Таким образом структура параллельного алгоритма включает в себя следующие этапы:

1. Начальное распределение вершин по процессорам. Начальное распределение может быть любым, например, в соответствии с порядковыми номерами вершин.
2. Определение числа доменов, которые должны быть сформированы на каждом из процессоров, и количества вершин в них.
3. Рекурсивная координатная бисекция вершин по процессорам. На каждом этапе сначала вычисляется координатная ось, вдоль которой область имеет наибольшую протяженность. Затем блок вершин бьется на две части перпендикулярно данной оси. Группа процессоров делится на две, далее каждая из групп делит свой блок вершин аналогичным образом. Для разделения блока вершин используется параллельная сортировка. После рекурсивной координатной бисекции вершин по процессорам на процессорах оказывается столько вершин, сколько должно быть в образуемых на них доменах.
4. Локальная рекурсивная координатная бисекция вершин по доменам. Дальнейшее разбиение на домены проводится локально на каждом процессоре.

Здесь приводится оценка параллельной сложности алгоритма: числа шагов, за которое можно выполнить данный алгоритм в предположении доступности неограниченного числа необходимых процессоров (функциональных устройств, вычислительных узлов, ядер и т.п.). Параллельная сложность алгоритма понимается как высота канонической ярусно-параллельной формы [1]. Необходимо указать, в терминах каких операций дается оценка. Необходимо описать сбалансированность параллельных шагов по числу и типу операций, что определяется шириной ярусов канонической ярусно-параллельной формы и составом операций на ярусах.

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

Для описания ресурса параллелизма алгоритма (ресурса параллелизма информационного графа) необходимо указать ключевые параллельные ветви в терминах конечного и массового параллелизма. Далеко не всегда ресурс параллелизма выражается просто, например, через координатный параллелизм или, что то же самое, через независимость итераций некоторых циклов (да-да-да, циклы - это понятие, возникающее лишь на этапе реализации, но здесь все так связано… В данном случае, координатный параллелизм означает, что информационно независимые вершины лежат на гиперплоскостях, перпендикулярных одной из координатных осей). С этой точки зрения, не менее важен и ресурс скошенного параллелизма. В отличие от координатного параллелизма, скошенный параллелизм намного сложнее использовать на практике, но знать о нем необходимо, поскольку иногда других вариантов и не остается: нужно оценить потенциал алгоритма, и лишь после этого, взвесив все альтернативы, принимать решение о конкретной параллельной реализации. Хорошей иллюстрацией может служить алгоритм, структура которого показана на рис.2: координатного параллелизма нет, но есть параллелизм скошенный, использование которого снижает сложность алгоритма с [math]n\times m[/math] в последовательном случае до [math](n+m-1)[/math] в параллельном варианте.

Рассмотрим алгоритмы, последовательная сложность которых уже оценивалась в п.1.6. Параллельная сложность алгоритма суммирования элементов вектора сдваиванием равна [math]\log_2n[/math], причем число операций на каждом ярусе убывает с [math]n/2[/math] до [math]1[/math]. Параллельная сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки - [math]\log_2n[/math]. Параллельная сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это [math]n[/math] шагов для вычислений квадратного корня, [math](n-1)[/math] шагов для операций деления и [math](n-1)[/math] шагов для операций умножения и сложения.

1.9 Входные и выходные данные алгоритма

Входные данные алгоритма: множество координат вершин [math]C = \{c_i\}[/math] графа [math]G = (V, E)[/math], где [math]c_i = (v_i^x, v_i^y, v_i^z)[/math] - координаты вершины [math]v_i \in V[/math]. Таким образом, множество координат вершин может быть передано как массив трёхмерных векторов. Число доменов [math]p[/math] передается в виде натурального числа, не может превышать число вершин.

Количество элементов массива входных данных может быть очень большим, например, алгоритм используется для нерегулярных сеток, содержащих [math]10^9[/math] и более вершин. В настоящее время такие сетки невозможно разместить в памяти одного процессора (на гексаэдральную сетку, состоящую из [math]1.2 \cdot 10^8[/math] ячеек, требуется порядка 200 ГБ), поэтому для декомпозиции нужен параллельный алгоритм.

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

Выходные данные алгоритма: разбиение множества вершин [math]V[/math] на заданное число [math]p[/math] доменов [math]S_j[/math]. Эти данные могут быть записаны как массив четырехмерных векторов координат вершин, четвертой координатой будет номер домена, которому принадлежит вершина.

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов является // TODO.

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

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

Так как основной операцией алгоритма является сравнение целых чисел при сортировке, наблюдается дисбаланс арифметических операций по отношению к операциям обращения к памяти. Алгоритм обладает сбалансированностью операций между параллельными ветвями алгоритма, что следует из условия [math](1)[/math] - количество вершин на разных процессах примерно одинаково.

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

"Длинные" дуги в информационном графе [1] говорят о потенциальных сложностях с размещением данных в иерархии памяти компьютера на этапе выполнения программы. С одной стороны, длина дуги зависит от выбора конкретной системы координат, в которой расположены вершины графа, а потому в другой системе координат они попросту могут исчезнуть (но не появится ли одновременно других длинных дуг?). А с другой стороны, вне зависимости от системы координат их присутствие может быть сигналом о необходимости длительного хранения данных на определенном уровне иерархии, что накладывает дополнительные ограничения на эффективность реализации алгоритма. Одной из причин возникновения длинных дуг являются рассылки скалярных величин по всем итерациям какого-либо цикла: в таком виде длинные дуги не вызывают каких-либо серьезных проблем на практике.

Для проектирования специализированных процессоров или реализации алгоритма на ПЛИС представляют интерес компактные укладки информационного графа [1], которые также имеет смысл привести в данном разделе.

2 Программная реализация алгоритма

2.4 Масштабируемость алгоритма и его реализации

Раздел будет добавлен позднее

2.7 Существующие реализации алгоритма

Алгоритм, как и многие другие алгоритмы разбиения графов реализованы в следующих последовательных пакетах декомпозиции графов: METIS, JOSTLE, SCOTCH, CHACO и PARTY. К параллельным пакетам относятся PARMETIS (параллельная версия пакета METIS), JOSTLE, PT-SCOTCH (параллельная версия пакета SCOTCH) и ZOLTAN.

ParMETIS и METIS распространяются свободно URL: http://glaros.dtc.umn.edu/gkhome/software .

ZOLTAN - свободная лицензия, URL: http://www.cs.sandia.gov/~web1400/1400_download.html .

Также алгоритм реализован в параллельном пакете GridSpiderPar, URL: http://lira.imamod.ru/FondProgramm/Decomposition/

3 Литература

[1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.

[2] Якобовский М.В. Введение в параллельные методы решения задач: Учебное пособие / Предисл.: В. А. Садовничий. – М.: Издательство Московского университета, 2012. – 328 с., илл. – (Серия «Суперкомпьютерное образование»), ISBN 978-5-211-06382-2 URL: http://lira.imamod.ru/ITTPMOPS/

[3] Головченко Е.Н. Параллельный пакет декомпозиции больших сеток // Математическое моделирование. 2011. Т. 23. № 10. 3-18. URL: http://cyberleninka.ru/article/n/razbienie-bolshih-setok

[4] Головченко Е.Н. Декомпозиция расчетных сеток для решения задач механики сплошных сред на высокопроиводительных вычислительных системах. Дисс ... канд. ф.-м.наук. Институт прикладной математики им. М. В. Келдыша РАН, г. Москва URL: http://keldysh.ru/council/3/D00202403/golovchenko_diss.pdf