[Из песочницы] Метод рекурсивной координатной бисекции для декомпозиции расчетных сеток08.02.2017 19:03

Введение

1. Повысить скорость работы программы.
2. Работать с сетками такого размера, который не помещается в оперативной памяти одного процессора.

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

Примечание: поскольку в алгоритме активно используется параллельная сортировка, то для понимания статьи настоятельно рекомендуется знать что такое «Сеть обменной сортировки со слиянием Бэтчера».

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

Сетки мы ограничим лишь регулярными двумерными, т.е. такими, в которых узел с индексами (i, j) соединен с соседними существующими по i, j узлами: (i-1, j), (i+1, j), (i, j-1), (i, j+1). Как нетрудно заметить, на верхней картинке сетка не регулярна, примеры регулярных сеток изображены ниже:

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

На вход программе поступают 3 аргумента:

• n1, n2 — размеры двумерной сетки.
• k — число доменов (частей) на которое требуется разбить сетку.

Один узел сетки мы будем хранить в следующей структуре:

struct Point{
    float coord[2]; // 0 - x, 1 - y
    int index;
};

Сама сетка хранится в массиве Point P[n1*n2], где узел с индексами (i, j) находится в P[i*n2+j]. В результате работы программы число вершин в доменах должно быть равно с точностью до одной вершины величине (n1*n2/k).

Алгоритм решения

1. Сортировка массива узлов вдоль одной из координатных осей (в двумерном случае вдоль x или y).
2. Разбиение отсортированного массива на 2 части.
3. Рекурсивный запуск процедуры для полученных подмассивов.

1. Рандомно.

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

2. Из геометрических соображений.

   Минимизируется длина проводимого разреза (в общем случае разрез — это гиперплоскость, поэтому в многомерном случае корректнее говорить о площади). Для этого нужно выбрать точки с минимальными и максимальными значениями координат после чего измерить длину на основании какой-либо метрики. Затем достаточно сравнить полученные длины и выбрать соответствующую ось. Конечно, здесь легко подобрать контрпример:

   

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

3. Минимизировать число разрезанных ребер.

   Разрезанное ребро — это ребро, соединяющее вершины из разных доменов. Таким образом, если минимизировать это число, то домены получаются как можно более «автономными», поэтому можно говорить о высоком качестве разбиения. Так, на картинке выше этот алгоритм выберет горизонтальную ось. Платить за это приходится скоростью, ведь нужно поочередно отсортировать массив вдоль каждой из осей, а затем посчитать число разрезанных ребер.

int k1 = (k + 1) / 2;
int k2 = k - k1;
int m1 = m * (k1 / (double)k);
int m2 = m - m1;

Примечание: некоторые могут спросить, зачем в выражении для m1 нужна вещественная арифметика, ведь без ее использования получится тот же результат? Все дело в переполнении, например, при разбиении сетки 10000×10000 на 100 доменов m * k1 = 1010, что выходит за границы диапазона int. Будьте внимательны!

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

void local_decompose(Point *a, int *domains, int domain_start, int k, int array_start, int n);

• a — исходная сетка.
• domains — выходной массив, в котором хранятся номера доменов для узлов сетки (имеет такую же длину что и массив a).
• domain_start — начальный индекс доступных номеров доменов.
• k — число доступных элементов.
• array_start — начальный элемент в массиве.
• n — число элементов в массиве.

void local_decompose(Point *a, int *domains, int domain_start, int k, int array_start, int n)
{
    if(k == 1){
        for(int i = 0; i < n; i++)
            domains[array_start + i] = domain_start;
        return;
    }

    axis = !axis;
    qsort(a + array_start, n, sizeof(Point), compare_points);

int compare_points(const void *a, const void *b)
{
    if(((Point*)a)->coord[axis] < ((Point*)b)->coord[axis]){
        return -1;
    }else if(((Point*)a)->coord[axis] > ((Point*)b)->coord[axis]){
        return 1;
    }else{
        return 0
    }
}

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

    local_decompose(a, domains, domain_start, k1, array_start, n1);
    local_decompose(a, domains, domain_start + k1, k2, array_start + n1, n2);
}

Параллельный алгоритм

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

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

Сначала разберемся со второй проблемой. Решение тривиально — нужно ввести фиктивные элементы, о которых вскользь упоминалось ранее. Вот тут нам и пригодится поле index в структуре Point. Делаем его равным -1 — и перед нами фиктивный элемент! Вроде все отлично, но этим решением мы существенно усугубили первую проблему. В общем случае на этапе разбиения массива на 2 части становится невозможно определить не только сам разбивающий элемент, но даже процессор, на котором он находится. Поясню это на примере: пусть у нас есть 9 элементов (сетка 3×3), которые нужно разбить на 3 домена на 4 процессорах, т.е. n = 9, k = 3, p = 4. Тогда после сортировки фиктивные элементы могут оказаться в любом месте (на картинке зеленым обозначены узлы сетки, а красным — фиктивные элементы):

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

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

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

int pc = n1 / elem_per_proc;

int middle = n1 % elem_per_proc;

Теперь процессоры с номерами меньше чем pc будут дальше работать с первой частью исходного массива, а остальные, параллельно, со второй частью. При этом в рамках одной группы число элементов на всех процессорах одинаковое (хотя у второй группы оно может быть другим), что позволяет дальше использовать сортировку Бетчера. Не стоит забывать, что возможен случай pc = 0 — тогда достаточно пересылать элементы в «другую сторону», т.е. на процессоры с большими номерами.

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

С учетом всего вышесказанного общая схема функции будет выглядеть так:

void decompose(Point **array, int **domains, int domain_start, int k, int n,
        int *elem_per_proc, MPI_Comm communicator)
{
    // инициализация переменных
    int rank, proc_count;
    MPI_Comm_rank(communicator, &rank);
    MPI_Comm_size(communicator, &proc_count);
    // базис рекурсии
    if(proc_count == 1){
        // если закончились процессоры, то запустить последовательный алгоритм
        local_decompose(...);
        return;
    }
    if(k == 1){
        // а если закончились домены, то присвоить оставшийся домен всем доступным элементам
        return;
    }

    // разделение массива на 2 части
    int k1 = (k + 1) / 2;
    int k2 = k - k1;
    int n1 = n * (k1 / (double)k);
    int n2 = n - n1;
    // нахождение разбивающего элемента
    int pc = n1 / (*elem_per_proc);
    int middle = n1 % (*elem_per_proc);
    // разделение процессоров на 2 группы, которые будут параллельно
    // работать с новыми подмассивами
    int color;
    if(pc == 0)
        color = rank > pc ? 0 : 1;
    else
        color = rank >= pc ? 0 : 1;
    MPI_Comm new_comm;
    MPI_Comm_split(communicator, color, rank, &new_comm);

    // выбор оси для сортировки и сама параллельная сортировка
    axis = ...
    sort_array(*array, *elem_per_proc, communicator);

    if(pc == 0){
        // перераспределение элементов на процессоры с большими номерами
        // и соответствующий рекурсивный вызов
        return;
    }
    // перераспределение элементов на процессоры с меньшими номерами
    // и немного другой рекурсивный вызов
    if(rank < pc)
        decompose(array, domains, domain_start, k1, n1, elem_per_proc, new_comm);
    else
        decompose(array, domains, domain_start + k1, k2, n2, elem_per_proc, new_comm);
}

Заключение

Комментарии (0)