Метод рекурсивну координатної бісекції для декомпозиції розрахункових сіток



Введення
Розрахункові сітки широко застосовуються при вирішенні чисельних задач за допомогою методів кінцевих різниць. Якість побудови такої сітки в значній мірі визначає успіх у рішенні, тому іноді сітки досягають величезних розмірів. У цьому випадку на допомогу приходять багатопроцесорні системи, адже вони дозволяють відразу вирішити 2 завдання:

  1. Підвищити швидкість роботи програми.
  2. Працювати з сітками такого розміру, який не поміщається в оперативній пам'яті одного процесора.
При такому підході сітка, що покриває розрахункову область, розбивається на безліч доменів, кожен з яких обробляється окремим процесором. Основна проблема тут полягає в «чесності» розбиття: потрібно вибрати таку декомпозицію, при якій обчислювальне навантаження розподілено рівномірно між процесорами, а накладні витрати, викликані дублюванням обчислень і необхідністю передачі даних між процесорами, малі.

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

Примітка: оскільки в алгоритмі активно використовується паралельна сортування, то для розуміння статті настійно рекомендується знати що таке "Мережа обмінної сортування злиттям Бэтчера".

Постановка завдання
Для простоти будемо вважати, що процесори у нас одноядерні. Ще раз, щоб не було плутанини:
1 процесор = 1 ядро. Обчислювальна система з розподіленою пам'яттю, тому будемо використовувати технологію MPI. На практиці в таких системах багатоядерні процесори, а значить для максимально ефективного використання обчислювальних ресурсів слід також використовувати треди.

Сітки ми обмежимо лише регулярними двовимірними, тобто такими, у яких вузол з індексами (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;
};

Крім координат точки тут ще є індекс вузла, який має 2 призначення. По-перше, його достатньо для відновлення топології сітки. По-друге, за допомогою нього можна відрізняти фіктивні елементи від нормальних (наприклад, у фіктивних елементів значення цього поля встановити рівним -1). Про те, що це за фіктивні елементи і навіщо вони потрібні, буде детально розказано далі.

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

Алгоритм вирішення
Процедура рекурсивну координатної бісекції складається з 3 кроків:

  1. Сортування масиву вузлів вздовж однієї з координатних осей (у двовимірному випадку уздовж x або y).
  2. Розбиття відсортованого масиву на 2 частини.
  3. Рекурсивний запуск процедури для отриманих подмассивов.
Базисом рекурсії тут є випадок k = 1 (залишився 1 домен). Відбувається приблизно наступне:

А може і так:


Звідки взялася неоднозначність? Неважко помітити, що я ніяк не обмовив критерій вибору осі, саме звідси і з'явилися різні варіанти (на 1 крок у першому випадку сортування проходить по осі x, а в другому по y). Так як вибрати вісь?

  1. Рандомно.

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

  2. З геометричних міркувань.

    Мінімізується довжина проведеного розрізу (в загальному випадку розріз — це гіперплощина, тому в багатовимірному випадку коректніше говорити про площі). Для цього потрібно вибрати точки з мінімальними і максимальними значеннями координат після чого виміряти довжину на підставі якої-небудь метрики. Потім достатньо порівняти отримані довжини і вибрати відповідну вісь. Звичайно, тут легко підібрати контрприклад:



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

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

    Розрізане ребро — це ребро, що з'єднує вершини з різних доменів. Таким чином, якщо мінімізувати це число, то домени виходять як можна більш «автономними», тому можна говорити про високий як розбиття. Так, на картинці вище цей алгоритм вибере горизонтальну вісь. Платити за це доводиться швидкістю, адже потрібно по черзі відсортувати масив вздовж кожної з осей, а потім порахувати число розрізаних ребер.
Ще варто обумовити як розбивати масив. Нехай m = n1*n2 (загальне число вузлів і довжина масиву P), k — число доменів (як і раніше). Тоді ми намагаємося розділити число доменів навпіл, а потім ділимо масив в такому ж співвідношенні. У вигляді формул:

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

На прикладі: розбити 9 елементів на 3 домену (m = 9, k = 3). Тоді спочатку масив розділиться в співвідношенні 6 до 3, а потім масив з 6 елементів розіб'ється навпіл. В результаті ми отримаємо 3 домену за 3 елемента в кожному. Те що треба.

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

Послідовний алгоритм
Розібравшись з теорією, можна переходити до реалізації. Для сортування будемо використовувати функцію qsort() із стандартної бібліотеки мови Сі. Наша функція буде приймати 6 аргументів:

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;
}

Якщо у нас залишився всього 1 домен, то ми поміщаємо всі доступні вузли сітки у нього. Потім ми вибираємо вісь і сортуємо сітку вздовж неї:

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

Вісь можна вибрати один із запропонованих вище способів. У цьому прикладі, як говорилося раніше, вона просто змінюється на протилежну. Тут axis — це глобальна змінна, яка використовується у функції порівняння елементів 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);
}

От і все. Повністю вихідний код функції local_decompose() доступний на гітхабі.

Паралельний алгоритм
Незважаючи на те, що алгоритм залишився тим же, його паралельна реалізація куди складніше. Я пов'язую це з двома основними причинами:

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;

Коли розбиває елемент визначено, що містить його процесор починає процедуру перерозподілу елементів: всі елементи до розбиває він рівномірно розсилає процесорам «до» нього (тобто з номерами, меншими ніж у нього), при необхідності додаючи фіктивні. У себе відправлені елементи він замінює на фіктивні. В картинках для n = 15, k = 5, p = 4:




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

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

З урахуванням всього вищесказаного загальна схема функції буде виглядати так:

Загальна схема функції decompose
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 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.