Участник:Fokina/Рекурсивная координатная бисекция

Рекурсивная координатная бисекция
Последовательный алгоритм
Последовательная сложность	[math]O((k-1)*n*log(n))[/math]
Объём входных данных	[math]n*m[/math]
Объём выходных данных	[math]n*m[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(log_2 k)[/math]
Ширина ярусно-параллельной формы	[math]O(k/2)[/math]

Автор: Фокина Н.Ю.

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

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

Метод рекурсивной координатной бисекции был предложен Бергером (M. Berger) и Бохари (S. Bokhari) в 1987 году для решения задачи статической баланисировки загрузки. В основе метода лежит идея равномерного разбиения точек в пространстве с помощью гиперплоскостей.

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

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

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

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

В начале работы имеется домен -- множество вершин графа [math]G[/math], заданных в виде векторов координат [math](v_1,..,v_m)[/math], где [math]m[/math] -- размерность пространства. На каждой итерации осуществляется сортировка вершин, входящих в домен, вдоль оси, по которой он имеет наибольшую протяженность. Затем домен разделяется на поддомены путем разбиения: первые [math]n/2[/math] или [math]n*k/2[/math] из отсортированного ранее домена относят к первому поддомену, оставшиеся -- к второму. Полученные поддомены также являются доменами и передаются в качестве входных данных на следующую итерацию до тек пор, пока общее количество поддоменов не достигнет необходимого числа [math]k[/math].

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

Основное время работы алгоритма приходится на сортировку данных на каждой итерации. В связи с этим целесообразно использовать алгоритмы сортировки, имеющие наименьшую среднюю сложность и показавшие высокую эффективность на практике. В последовательном случае такой сортировкой может являться комбинация из пирамидальной сортировки и сортировки слиянием, имеющая сложность [math]O(n*\log(n))[/math]. Параллельная версия данной сортировки также эффективна и имеет сложность [math]O(\frac{n}{p} \big\lbrack \log(\frac{n}{p})+\frac{\lceil \log (p) \rceil^2}{2} \big\rbrack)[/math]. Еще больший прирост производительности может дать битоническая сортировка, сложность которой оценивается как [math]O(\log^2(n))[/math].

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

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

На каждом этапе рекурсивной координатной бисекции происходит следующее:

1. Определяются границы наименьшего параллелепипеда, охватывающего все вершины, которые делятся на данном этапе. Вычисляются максимумы по каждой из координат. Определяется координата [math]j[/math], вдоль которой параллелепипед имеет наибольшую протяженность.
2. Выполняется сортировка вершин по координате [math]j[/math].
3. Создается два новых поддомена. Первая половина вершин из отсортированного домена помещается в первый поддомен, все последующие -- в второй.

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

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

Реализация метода рекурсивной бисекции на языке C в общем виде изображена на следующем рисунке:

void
recursive_bisect(struct Vertex *graph, int n, int k) {
    int                             k1, k2;
    int                             n1, n2;
    struct Vertex *                 subgraph1;
    struct Vertex *                 subgraph2;

    k1 = (k + 1) / 2;
    k2 = k - k1;
    n1 = n * k1 / k;
    n2 = n - n1;

    bisect(graph, n, subgraph1, n1, subgraph2, n2);
    
    if (k1 > 1) recursive_bisect(subgraph1, n1, k1);
    if (k2 > 1) recursive_bisect(subgraph2, n2, k2);

    return;
}

В качестве входных параметров функция принимает массив вершин, содержащихся в текущем домене (graph), числа n -- количество вершин, k -- количество поддоменов, на который требуется разбить данный домен.

В теле функции вычисляется отношение между количеством поддоменов, на которые будут впоследствии разбиты полученные на текущей итерации поддомены (k1 и k2), и размеры поддоменов, которые будут получены на данной итерации (n1 и n2). С помощью вызова функции bisect осуществляется непосредственно разбиение на два поддомена. Затем функция recursive_bisect рекурсивно вызывается на полученных поддоменах до тех пор, пока остается необходимость в дальнейшем их разбиении.

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

void
coord_bisect(struct Vertex *graph, int n, struct Vertex *subgraph1, int n1, struct Vertex *subgraph2, int n2) {
    int                             axis;

    /* Calculate direction in which the initial graph is most elongated. */
    axis = calculate_max_distance(graph, n);

    /* Sort vertices along the selected axis. */
    sort(graph, n, axis);

    /* Create subgraphs. */
    subgraph1 = create_subgraph(graph, 0, n1);
    subgraph2 = create_subgraph(graph, n1, n);

    return;
}

Входными параметрами функции являются массив вершин в домене, его размер, а также размеры полученных в результате разбиения поддоменов. Выходными параметрами являются массивы вершин, входящих в полученные поддомены. Следует отметить, что на вход достаточно передавать размер лишь одного из поддоменов, например, первого, поскольку размер второго явно из него следует, т.к. [math]n1 + n2 == n[/math]. В данном примере оба размера передаются для наглядности и единообразия.

В процессе работы функции с помощью вызова функции calculate_max_distance выбирается координатная ось, вдоль которой сетка имеет наибольшую протяженность. Затем вершины сортируются вдоль этой оси, и создаются два поддомена: в первый входят вершины с индексами [0..n1) в отсортированном массиве, в второй -- [n1..n), где n1 равно [math]\frac{n*\frac{(k+1)}{2}}{2}[/math].

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

В течение каждого этапа выполняются следующие действия:

1. Определение диаметра графа по каждой из координат -- [math]m*O(n)[/math], где [math]m[/math] -- размерность пространства, [math]n[/math] -- количество вершин. Поскольку величина [math]m[/math], как правило, очень мала, ей можно пренебречь, оставив в качестве оценки сложности данного действия [math]O(n)[/math].
2. Поиск максимального расстояния -- осуществляется за [math]O(m)[/math], также пренебрежимо мало.
3. Сортировка. Наиболее эффективным алгоритмом последовательной сортировки является комбинация сортировки слиянием и пирамидальной сортировки . Сложность как пирамидальной сортировки, так и сортировки слиянием оценивается как [math]O(nlog_2(n))[/math].
4. Разбиение домена на поддомены осуществляется за [math]O(1)[/math].

Исходя из этого, общая сложность каждой итерации оценивается сверху как [math]O(nlog_2(n))[/math], а их количество составляет [math](k - 1)[/math]. Таким образом, последовательная сложность метода -- [math]O((k-1) *n*log(n))[/math].

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

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

Рис. 1. Общая схема рекурсивной координатной бисекции

На рисунке 2 изображен информационный граф метода при распределении домена из 4 вершин, заданных в двумерном пространстве, на 2 поддомена. Входными данными являются вектора A1-A4, которые могут содержаться либо в исходном домене, либо в поддомене, полученном на предыдущей итерации. Выходные данные -- непересекающиеся множества векторов B1 и B2. В них содержатся вектора A1-A4. Вершины с меткой CD обозначают операцию вычисления диаметра домена (т.е. наибольшего расстояния между содержащимися в домене вершинами) по координатам [math]x[/math] и [math]y[/math] соответственно. Вершина MAX обозначает нахождение максимума, SORT -- сортировку, SPLIT -- разбиение домена на два поддомена равного размера.

Рис. 2. Информационный граф для одной итерации

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

При разбиении графа на [math]k[/math] подграфов глубина рекурсии составляет [math]\lceil log_2 k \rceil[/math], а на каждом [math]i[/math] шаге решается до [math]2^i[/math] подзадач. Поскольку подзадачи, возникающие на каждом шаге рекурсии, не имеют информационных зависимостей друг между другом, их решение можно распараллелить, получив таким образом некоторый прирост производительности. Исходя из того, что сложность этих подзадач примерно одинакова, можно считать, что их решение занимает примерно одинаковое время [math]t[/math] и не зависит от яруса, на котором решается конкретная подзадача. При разбиении графа на [math]k[/math] подграфов требуется решить [math]k-1[/math] подзадач. Время их последовательного выполнения составляет [math]T_1 = t*(k-1)[/math]. Пусть доступно бесконечное число одинаковых процессоров. При распараллеливании затраченное время сокращается до [math]T_\infty(k) = t * \lceil log_2 k \rceil[/math]. Таким образом, предельное ускорение от распараллеливания рекурсивной бисекции составляет [math]S(k) = \frac{T_1(k)}{T_\infty(k)} = \frac{k-1}{\lceil log_2 k \rceil}[/math].

Если количество доступных процессоров ограничено некоторым числом [math]p[/math], то имеет место следующее соотношение: [math]S(k, p) = \frac{p * \lfloor log_2 k \rfloor}{2*p - 2 + \lfloor log_2 k \rfloor - \lceil log_2 p \rceil}[/math].

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

Метод рекурсивной координатной бисекции работает над графами, часто -- сеточными. Поскольку метод не учитывает связи между вершинами, входные данные удобно задавать в виде массива координат в [math]n[/math]-мерном прострастве, т.е. векторов размерности [math]n[/math]. Как правило, в прикладных задачах [math]n[/math] составляет 2 или 3, т.е. заданы координаты на плоскости или в пространстве.

Входные данные алгоритма: Исходный граф [math]G[/math], натуральное число [math] k \in \N [/math] -- количество подграфов, на которое нужно осуществить разбиение.

Объем входных данных: [math]n*m[/math], где [math]n[/math] -- количество векторов, [math]m[/math] -- размерность пространства.

Выходные данные алгоритма: Множество [math]G_1 .. G_k[/math] подграфов исходного графа [math]G[/math].

Объем выходных данных: [math]n*m[/math], где [math]n[/math] -- количество векторов, [math]m[/math] -- размерность пространства.

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

Метод рекурсивной координатной бисекции является устойчивым и детерминированным. Его вычислительная мощность пропорциональна [math]log(n)[/math]. Вычислительная сложность определяется, главным образом, использованным алгоритмом сортировки и составляет [math]O((k-1)*n*log(n))[/math] в последовательном и [math]O(\log(k)*\log^2(n))[/math] в параллельном случае.

Можно также отметить следующие особенности метода:

• Использование параллельного алгоритма рекурсивной координатной бисекции целесообразно при разбиении на большое число подграфов.

• Наибольшее ускорение достигается при разбиении на число подграфов, являющееся степенью двойки.

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

• Небольшие накладные расходы, возникающие при обращении процессора к полной информации о разбиении.

• Отсутствие явного контроля над накладными расходами на передачу данных.

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

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

Масштабируемость алгоритма была исследована на примере реализации метода рекурсивной координатной бисекции, содержащейся в библиотеке ZOLTAN v3.83 [1] (файл example/C/simpleRCB.c). В части испытаний использованная программа была модифицирована таким образом, чтобы фиксировать количество поддоменов, на которое требуется разбить сетку. Для этого в функцию, непосредственно осуществляющую бисекцию Zoltan_LB_Partition(), был дополнительно передан параметр NUM_GLOBAL_PARTS (часть контекста struct Zoltan_Struct *zz).

В качестве тестового стенда было использовано ЭВМ BlueGene/P ВМК МГУ. Первый этап тестирования выполнялся в соответствии с методикой испытаний:

• фиксировано значение параметра [math]k[/math] -- итоговое количество поддоменов (16);
• в качестве тестового набора входных данных была сгенерирована двумерная регулярная сетка размера [math]N*N[/math], где [math]N[/math] -- 1024, 2048, 3072, 4096, 5120.
• для каждого размера сетки выполнена серия запусков на различном количестве процессоров: 1, 2, 4, [8..128] c шагом 8 (8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 128).

Тестовая программа и сама библиотека ZOLTAN были скомпилированы с помощью mpicc 1.1 (gcc 4.1.2), параметры компиляции стандартные, приведены в Makefile/make everything. Для запуска использовался скрипт, одновременно ставящий в очередь все задачи с одинаковым размером входной сетки. Пример запуска одной задачи: mpisubmit.bg -n $3 --stdout result_$2_$1_$3.txt -w 00:15:00 simpleRCB.exe -- mesh$1.txt $2 $1; $1 -- размер сетки, $2 -- количество поддоменов, $3 -- количество процессоров.

Полученные результаты производительности показаны на рис. 3 и 4. На рис. 4 для наглядности использована логарифмическая шкала для количества процессоров, а также исключены результаты для запусков при количестве процессоров, не являющемся степенью 2. Соответствующие им результаты эффективности распараллеливания приведены на рис. 5 и 6.

Рис. 3. Производительность алгоритма в зависимости от числа процессора и размера сетки

Рис. 4. Производительность алгоритма в зависимости от числа процессора и размера сетки (упрощенно)

Рис. 5. Эффективность распараллеливания алгоритма в зависимости от числа процессора и размера сетки

Рис. 6. Эффективность распараллеливания алгоритма в зависимости от числа процессора и размера сетки (упрощенно)

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

Минимальное значение эффективности распараллеливания составило 0.0169681072 (сетка 1024*1024, 128 процессоров), максимальное -- 0.7165250222 (сетка 4096*4096, 2 процессора).

На рис. 3 было обнаружено небольшое увеличение производительности алгоритма при количестве процессоров, равном количеству поддоменов, на которое требовалось разбить сетку (16). Это обусловлено тем, что в случае совпадения количества поддоменов и процессоров, количество вершин, приходящихся на каждый из процессоров в течение обработки, оптимально. В связи с этим было проведено дополнительное исследование производительности алгоритма в таком режиме, т.е. при всех запусках итоговое количество поддоменов совпадало с количеством используемых процессоров. Такое использование характерно для решения задачи балансировки нагрузки, а также предварительного разбиения вершин сетки по процессорам, для чего обычно и используется рекурсивная координатная бисекция. Результат тестирования приведен на рис. 7.

Рис. 7. Производительность алгоритма в зависимости от числа процессора и размера сетки, количество поддоменов совпадает с количеством процессоров

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

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

Наиболее широко в настоящее время применяется параллельный алгоритм рекурсивной координатной бисекции, реализованный в пакете ZOLTAN (Zoltan2) [2].

Метод рекурсивной координатной бисекции сеток реализован в пакете GridSpiderPar [3]. Использованный в нем алгоритм отличается от аналогичного алгоритма в пакете ZOLTAN тем, что в нем секущая плоскость (медиана) при необходимости разрезается по нескольким координатам. Это позволяет обрабатывать ситуации наличия на одной плоскости множества узлов с одинаковым значением координаты, что характерно для регулярных сеток. В пакете ZOLTAN вершины из медианы распределяются по областям произвольным образом, что увеличивает число разрезанных ребер.

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