=Paper=
{{Paper
|id=Vol-1631/105-112
|storemode=property
|title=Research and development of Johnson's algorithm parallel schemes in GPGPU technology
|pdfUrl=https://ceur-ws.org/Vol-1631/105-112.pdf
|volume=Vol-1631
|authors=Sergiy Pogorilyy,Maxim Slynko
|dblpUrl=https://dblp.org/rec/conf/ukrprog/PogorilyyS16
}}
==Research and development of Johnson's algorithm parallel schemes in GPGPU technology==
https://ceur-ws.org/Vol-1631/105-112.pdf
Proceedings of the 10th International Conference of Programming UkrPROG’2016 (Kyiv, Ukraine)
УДК 004.3
СТВОРЕННЯ І ДОСЛІДЖЕННЯ ПАРАЛЕЛЬНИХ СХЕМ
АЛГОРИТМУ ДЖОНСОНА В ТЕХНОЛОГІЇ GPGPU
С.Д. Погорілий, М.С. Слинько
Запропоновано застосування алгоритму Джонсона для знаходження найкоротших шляхів між усіма парами вершин зваженого
орієнтованого графа. Виконано його формалізацію у термінах модифікованих систем алгоритмічних алгебр Глушкова. Обґрунтова-
но доцільність використання технології GPGPU для пришвидшення роботи алгоритму. Отримано низку схем паралельної версії ал-
горитму, оптимізовану для використання в технології GPGPU. Запропоновано підходи до реалізації отриманих схем з використан-
ням архітектури обчислень NVIDIA CUDA. Виконано експериментальне дослідження підвищення продуктивності при проведенні
обчислень на відеоадаптері.
Ключові слова: NVidia CUDA, GPGPU, SSSP, APSP, Thrust, САА схема, алгоритм Джонсона.
Предложено применение алгоритма Джонсона для нахождения кратчайших путей между всеми парами вершин взвешенного ори-
ентированного графа. Выполнена его формализация в терминах модифицированных систем алгоритмических алгебр Глушкова.
Обоснована целесообразность использования технологии GPGPU для ускорения работы алгоритма. Получен ряд схем параллель-
ной версии алгоритма, оптимизированной под использование в технологии GPGPU. Предложено подходы к реализации получен-
ных схем с использованием архитектуры вычислений NVIDIA CUDA. Выполнено экспериментальное исследование повышения
производительности при проведении вычислений на видеоадаптере.
Ключевые слова: NVidia CUDA, GPGPU, SSSP, APSP, Thrust, САА схема, алгоритм Джонсона.
Johnson’s all pairs shortest path algorithm application in an edge weighted, directed graph is considered. Its formalization in terms of
Glushkov’s modified systems of algorithmic algebras was made. The expediency of using GPGPU technology to accelerate the algorithm is
proved. A number of schemas of parallel algorithm optimized for using in GPGPU were obtained. Suggested approach to the implementation
of the schemes obtained using computing architecture NVIDIA CUDA. An experimental study of improved performance by using GPU for
computations was made.
Key words: NVidia CUDA, GPGPU, SSSP, APSP, Thrust, SAA scheme, Johnson’s algorithm.
Вступ
Багато застосувань теорії графів існують у галузі мережевого аналізу. Однією із найважливіших і найак-
туальніших у сьогоденні областей застосування задачі пошуку найкоротших шляхів у ґрафах є маршрутизація
пакетів у комп’ютерних мережах [1]. Теорія графів використовується в хімії та молекулярній біології для задач
моделювання. Графи – зручна модель для розв’язання таких задач, як дизайн VLSI мікросхем, філогенетика,
видобуток даних, біоінформатика, аналіз мереж тощо.
В роботі розглядається проблема знаходження найкоротших шляхів між усіма парами вершин зваженого
орієнтованого графа (APSP), яка потребує великих обчислювальних потужностей. Наприклад, алгоритм
Флойда – Уоршелла (1962) має складність, O(| V 2 |) , де | V | – кількість вершин у графі. Різні алгоритми пропо-
нують різні підходи до вирішення задачі: наприклад, алгоритм Данцига полягає у послідовному обчисленні за
допомогою рекурентної процедури підматриць найкоротших шляхів зростаючої розмірності. Втім, найбільш
очевидним шляхом є ітеративне застосування алгоритму знаходження всіх найкоротших шляхів від однієї вер-
шини (SSSP) для кожної з вершин графу. В роботі розглядається алгоритм Джонсона ([2]), який полягає у ітера-
тивному застосуванні алгоритму Дейкстри до кожної з вершин графу. Нотацію алгоритмів подано у системі
алгоритмічних алгебр Глушкова (САА).
Архітектури сучасних GPU від NVIDIA
Графічні адаптери спочатку використовувалися лише для обробки зображень. Основною особливістю
сучасних графічних адаптерів є наявність набору потокових мультипроцесорів (Streaming multiprocessor, SM),
які є універсальними пристроями обробки даних, що уможливлює широке використання GPU для обчислення
широкого кола задач (технологія GPGPU – General Purpose GPU).
CUDA (Computed Unified Device Architecture) – це архітектура паралельних обчислень, розроблена ком-
панією NVIDIA для спрощення програмування GPGPU за рахунок використання високорівневих API. В моделі
програмування CUDA CPU визначають як «хост», а GPU – як пристрій (обчислювальний) [3]. Фактично, при-
стрій CUDA є багатоядерним процесором, що використовується для обчислення задач, та який може бути заді-
яним лише з CPU. За таксономією Флінна, принцип роботи CUDA пристроїв є близьким до принципу SIMD
(Singe Instruction – Multiple Data), з деякими уточненнями: кожен потік не лише виконує одну й ту ж саму ін-
струкцію над своєю підмножиною даних, але й може мати власний шлях виконання в залежності від заданих
умов. Розробники CUDA називають такий підхід SIMT (Single Instruction, Multiple Thread), тобто одна інструк-
ція виконується одночасно багатьма потоками, при чому поведінка кожного окремого потоку нічим не обмежу-
ється. Таким чином, основною програмною обчислювальною одиницею є потік. Блок – це набір одночасно
виконуваних потоків, що можуть взаємодіяти за допомогою спільної пам’яті та бар’єрної синхронізації. Сітка –
це набір блоків, що виконують одну і ту ж саму функцію-ядро [3, 4].
105
Proceedings of the 10th International Conference of Programming UkrPROG’2016 (Kyiv, Ukraine)
CUDA визначає розширення мови C/C++, яке дозволяє створювати код, що виконуватиметься на відеоа-
даптері, визначаючи С функції-ядра (kernels). Кожен потік можна однозначно визначити за допомогою ряду
контекстних змінних.
1. gridDim – векторна змінна, що зберігає розмірності сітки.
2. blockIdx – векторна змінна, що зберігає індекс (3-вимірний) блоку в сітці.
3. blockDim – векторна змінна, що зберігає розмірність блоку.
4. threadIdx – векторна змінна, що зберігає індекс (3-вимірний) потоку в блоці.
5. warpSize – змінна, що зберігає розмір warp групи [5].
Наприклад, наступна функція-ядро додає 2 числа:
__global__ void add(int* a, int* b, int* c){
*c = *a + *b;
}
int main(void){
... //ініціалізація змінних, виділення пам’яті на GPU
add<<<1, 1>>>(a, b, c); //запуск 1-потокової 1-блокової програми на GPU
return 0;
}
Система Алгоритмічних Алгебр
Для формалізованого представлення алгоритмів функціонування абстрактної моделі ЕОМ В.М. Глушков
запропонував математичний апарат систем алгоритмічних (мікропрограмних) алгебр (САА).
Фіксована САА являє собою двоосновну алгебраїчну систему, основами якої є множина операторів і
множина умов [6]. Операції САА поділяються на логічні і операторні. До логічних відносяться узагальнені
булеві операції і операція лівого множення оператора на умову (призначена для прогнозування обчислювально-
го процесу), а до операторних – основні конструкції структурного програмування (послідовне виконання, цик-
лічне виконання тощо).
Теорема Глушкова: Для довільного алгоритму існує (в загальному випадку не єдина) САА, в якій цей ал-
горитм може бути представлений регулярною схемою. Тобто, якщо визначити основи САА для конкретного
алгоритму, можна представити цей алгоритм у вигляді схеми і проводити подальші трансформації й оптимізації
вже не з алгоритмом, а з його САА схемою. Деякі основні операції алгоритмічної алгебри та їх аналоги у про-
цедурних мовах програмування наведено в табл. 1, 2:
Таблиця 1. Сигнатура операцій САА
Сигнатура операцій алгоритмічної алгебри
Відповідні оператори мови Паскаль
(САА)
Композиція: А * В (або A× B ) A;B;
α-диз’юнкція α ( A ∨ B ) if α then A else B;
α-ітерація α {A} while α do A;
Обернена α-ітерація {A}α repeat A until not α
Таблиця 2. Сигнатура операцій САА-М
Сигнатура паралельних
Позначення операції
операцій САА-М
Асинхронна диз’юнкція A || B Бінарна операція паралельного виконання операторів А і В
•
(або A ∨ B ) на підструктурах певної моделі
Фільтрація u Унарна операція, що породжує оператори-фільтри
Бінарна операція
Синхронна диз’юнкція A ∨ B синхронного застосування
операторів A i B
106
Proceedings of the 10th International Conference of Programming UkrPROG’2016 (Kyiv, Ukraine)
Алгоритм Джонсона. Ідея алгоритму полягає у багаторазовому застосуванні алгоритму Дейкстри до
кожної з вершин графу. Таким чином можна отримати найкоротші шляхи між будь-якими парами вершин. Не-
доліком алгоритму Дейкстри є неможливість його застосування для графів з від’ємними вагами ребер. Якщо ж
в графі є ребра з від’ємними вагами (але відсутні цикли негативних ваг), можна спробувати звести задачу до
випадку невід’ємних ваг, замінивши функцію ваги w на нову функцію ваги ŵ таку, щоб.
1. Найкоротші шляхи не змінилися: для будь-якої пари вершин u, v ∈ V найкоротший шлях з u в v з точ-
ки зору функції w є також найкоротшим шляхом з точки зору ŵ і навпаки.
2. Ваги всіх ребер з точки зору функції ŵ є невід’ємними [2].
Введемо наступні позначення:
G (V , E ) – вхідний граф;
V – початкова розмірність графа;
isEdge(i, j ) – предикат, значенням якого буде істина, якщо існує ребро з вершини i в j ;
d ij – довжина найкоротшого знайденого шляху з i в j ;
wij – довжина ребра з вершини i в j ;
{hi }, i = 0, V + 1 – множина найкоротших відстаней від доданої вершини до всіх вершин графа;
setWay(i, j , val ) – встановлення значення (val) шляху d ij ;
setEdge(i, j , val ) – встановлення значення (val) ребра з вершини i в j ;
add (Q, i ) – додавання елементу i в колекцію Q;
erase(Q, i ) – видалення елементу i з колекції Q;
eraseAll (Q, u ) – видалення всіх елементів колекції u з колекції Q;
size(Q) – отримання кількості елементів у колекції Q.
Тому першим кроком алгоритму є побудова графа G ′ = (V ′, E ′) , де
′
V = V ∪ {s} , а s – деяка нова вершина. При цьому мають бути створені ребра
{( s, v) : v ∈ V ∪ wsv = 0} . Оператор extendGraph розширює заданий граф додатковою вершиною (з порядковим
індексом V), додаючи ребра вагою 0 від доданої вершини до усіх існуючих:
extendGrap h = i d qu + wuv ( setWay(q,v ,d qu + wuv ) ∨ E ) ∨ E ) ; (3)
оператор RelaxEdges проводить релаксацію всіх ребер у графі:
RelaxEdges ( q ) = u Q(V);
INIT<<>>(Q, M, C, q);
while(Q.size() != 0)
{
RELAX_WITH_MASK<<>>(G, M, C);
threshold = ExtractMin(Q, C);
UPDATE_MASK<<>>(Q, C, M, threshold);
}
109
Proceedings of the 10th International Conference of Programming UkrPROG’2016 (Kyiv, Ukraine)
На кожній ітерації релаксація вихідних ребер відбувається не для всіх вершин, а з певною маскою:
RELAX_WITH_MASK(G, M, C)
{
tid = threadx.Id;
if(M[tid] == true)
{
<для всіх ребер [edgeId] що виходять з вершини tid виконати>
atomicMin(C[edgeId], G[edgeId] + C[tid]);
}
}
Як показано вище, вершини переносяться у масив встановлених, якщо відстань від початкової вершини
до них є меншою за поріг:
UPDATE_MASK(Q, C, M, threshold)
{
tid = thread.Id;
if(C[Q[tid]] <= threshold)
{
M[Q[tid]] = true;
Q[tid] = -1; // після виконання елементи ‘-1’ будуть видалені
}
}
Розробникам застосувань під GPU потрібно вирішувати не лише алгоритмічні задачі, але й задачі вирі-
шення конфліктів банків, оптимізації розгалужень тощо. При досить складних задачах кроки оптимізації є не-
очевидними та погіршують читабельність та зрозумілість коду. Розробники CUDA створили бібліотеку шабло-
нів на С++ для CUDA GPU застосувань, аналогічну до STL. Thrust вирішує низку проблем, пов’язаних з опти-
мізацією застосувань під конкретні архітектури, або ж використання різних підходів для пристроїв, що підтри-
мують різні версії CUDA. Розробники можуть повністю концентруватися на вирішенні високорівневих задач,
делегуючи Thrust вибір обчислювальної реалізації.
Основними контейнерами Thrust є вектори: host_vector та device_vector. Вміст host_vector зберігається в
загальній пам’яті комп’ютера, а вміст device_vector зберігається в пам’яті GPU. Контейнери реалізують низку
стандартних операцій та дають можливість глибокого розширення [10].
З формальної схеми видно, що у векторі Q зберігаються індекси невстановлених вершин. Алгоритм Дей-
кстри на кожній ітерації потребує знаходження індексу вершини, відстань від початкової вершини до якої є
найменшою. За допомогою Thrust задача вирішується наступним шляхом:
struct compare_unsettled_value
{
const int* data;
compare_unsettled_value(int* ptr) { data = ptr; }
__host__ __device__
bool operator()(int lhs, int rhs){
return this->data[lhs] < this->data[rhs];
}
};
ExtractMin(Q, C) //Q – device_vector
{
...
thrust::detail::normal_iterator> minIndex = thrust::min_element(
thrust::device, Q.begin(), Q.end(), compare_unsettled_value(C);
...
}
У вищенаведеному блоці коду використовується бібліотечний метод min_element, що виконує згортку
масиву (заданого межами Q.begin() – Q.end()) операцією знаходження мінімуму. Оскільки в масиві зберігаються
не значення відстані, а індекси, вводиться структура compare_unsettled_value, яка за заданими індексами прово-
дить порівняння відстаней відповідних вершин. Всі обчислення проводяться на GPU, про що свідчить явно
задана політика thrust::device [10].
110
Proceedings of the 10th International Conference of Programming UkrPROG’2016 (Kyiv, Ukraine)
Паралелізація внутрішнього циклу. Як було сказано вище, паралелізація внутрішнього циклу означає
одночасну релаксацію всіх ребер обраної вершини. Таким чином, оператор (13) залишається такою ж, як і в
послідовній версії ((8)), але змінюється логіка роботи з оператором (3): наразі кожен потік представляє собою
ребро, яке може бути релаксовано незалежно. Після паралельної обробки відбувається бар’єрна синхронізація
по V потокам, після якої обробка поточної вершини вважається закінченою і алгоритм обирає вершину для
наступної ітерації.
V
•
PSubD 2(u ) = ∨ ( RelaxEdge(q, u, v) × Bv (v)),
v =1
Dijkstra ( q ) = init ( q )×Q ≠∅ {ExtractMin (Q , u ) × PSubD 2(u )}. (16)
Паралелізація алгоритму Беллмана – Форда. Оскільки алгоритм Беллмана – Форда виконується в ал-
горитмі Джонсона лише раз, навіть при ідеальній паралельній версії великого виграшу отримати не вдасться.
Втім, алгоритм може опрацьовувати кожну вершину паралельно, тому оператор (4) у паралельному виконанні
має наступний вигляд:
V
•
RelaxEdges ( q ) = ∨ ( u