Запропоновано застосування алгоритму Джонсона для знаходження найкоротших шляхів між усіма парами вершин зваженого орієнтованого графа. Виконано його формалізацію у термінах модифікованих систем алгоритмічних алгебр Глушкова. Обґрунтовано доцільність використання технології 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.
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} Обернена α-ітерація {A}α Таблиця 2. Сигнатура операцій САА-М Відповідні оператори мови Паскаль
A;B; if α then A else B;
while α do A; repeat A until not α Синхронна диз’юнкція A ∨ B Позначення операції Бінарна операція паралельного виконання операторів А і В
на підструктурах певної моделі Унарна операція, що породжує оператори-фільтри
Бінарна операція синхронного застосування операторів A i B Алгоритм Джонсона. Ідея алгоритму полягає у багаторазовому застосуванні алгоритму Дейкстри до кожної з вершин графу. Таким чином можна отримати найкоротші шляхи між будь-якими парами вершин. Недоліком алгоритму Дейкстри є неможливість його застосування для графів з від’ємними вагами ребер. Якщо ж в графі є ребра з від’ємними вагами (але відсутні цикли негативних ваг), можна спробувати звести задачу до випадку невід’ємних ваг, замінивши функцію ваги w на нову функцію ваги wˆ таку, щоб.
1. Найкоротші шляхи не змінилися: для будь-якої пари вершин u, v ∈V найкоротший шлях з u в v з точки зору функції w є також найкоротшим шляхом з точки зору wˆ і навпаки.
2. Ваги всіх ребер з точки зору функції wˆ є невід’ємними [2]. Введемо наступні позначення:
G(V , E) – вхідний граф; V – початкова розмірність графа; isEdge(i, j) – предикат, значенням якого буде істина, якщо існує ребро з вершини i в j ; dij – довжина найкоротшого знайденого шляху з i в j ; wij – довжина ребра з вершини i в j ; {hi}, i = 0,V +1 – множина найкоротших відстаней від доданої вершини до всіх вершин графа; setWay(i, j, val) – встановлення значення (val) шляху dij ; 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. Тому першим кроком
алгоритму V ′ = V ∪{s} , а s – деяка нова вершина.
є При побудова цьому графа
G′ = (V ′, E′) , мають бути створені ребра {(s,v) : v ∈V ∪ wsv = 0} . Оператор extendGraph розширює заданий граф додатковою вершиною (з порядковим індексом V), додаючи ребра вагою 0 від доданої вершини до усіх існуючих:
extendGrap h=i<V {setEdge (V , i,0)} . Для формалізації алгоритмів введемо такі оператори: оператор init проводить ініціалізацію перед виконанням алгоритмів Дейкстри та Беллмана – Форда, встановлюючи найкоротші відстані від початкової вершини ( q ) до усіх інших вершин, окрім q , як ∞ : init(q)=i<V {i=q (setWay(q,i,0) ∨setWay(q,i,∞)) × (+ + i)} ; оператор RelaxEdge проводить релаксацію ребра з вершини u в v. Умовою релаксації є мінімум шляху від початкової вершини q до вершин, що розглядаються:
RelaxEdge (q,u,v)=isEdge(u,v)(duv >dqu +wuv (setWay(q,v ,d qu + wuv ) ∨ E) ∨ E) ; оператор RelaxEdges проводить релаксацію всіх ребер у графі:
RelaxEdges (q)=u<V {v<V {RelaxEdge (q,u,v ) × (+ + v)} × (+ + u)}.
Алгоритм Беллмана – Форда. Для коректного виконання алгоритму вхідний граф не може мати цикли негативної ваги. Для перевірки наявності такого циклу використовується алгоритм Беллмана – Форда для вершини s . Масив найкоротших шляхів від вершини s до усіх інших вершин, отриманий як результат виконання алгоритму Беллмана – Форда (який буде містити елементи, відмінні від 0 у випадку негативних ваг ребер), дозволяє перейти до вагової функції wˆ .
Згідно алгоритму Беллмана – Форда, після (V − 1) циклів релаксації усіх ребер встановлені шляхи duv є найкоротшими. Тому виконання умови релаксації після (V − 1) циклів релаксації усіх ребер означає наявність в графі цикла негативних ваг [2], що і перевіряє умова cnc:
u<V v<V cnc(q) = ∧ ∧ (isEdge(u, v) ∧ dqv < dqu + wuv ).
u=1 v=1
де
(
BellmanFor d (q) = init (q)×i<V −1 {Re laxEdges (q) × (+ + i)} × cnc(q).
Результатом виконання алгоритму Беллмана – Форда є булеве значення, що вказує на наявність циклів негативної ваги у графі. Також алгоритм визначає найкоротші шляхи від доданої вершини “q” до усіх існуючих вершин графа (задає значення {hi} ).
Якщо у графі немає циклів негативної ваги, можна трансформувати ваги його ребер:
UpdateWeights=i<V +1{ j<V +1 {isEdge(i, j) (wij = (wij + hi − h j ) ∨ E) × (+ + j)}× (+ + i)}.
Алгоритм Дейкстри. Алгоритм Дейкстри є по суті послідовним і знаходить найкоротші шляхи від заданої вершини до всіх інших вершин у графі за O(V logV + E) час. В алгоритмі Дейкстри виділяють 2 типи вершин: встановлені та невстановлені. Встановленими вершинами називають ті, для яких визначений шлях від початкової вершини і всі вихідні ребра яких релаксуються (або вже релаксовані). На початку алгоритму до масиву встановлених вершин додається початкова вершина, і релаксуються ребра, що з неї виходять. На наступній ітерації визначається вершина U, довжина до якої від початкової вершини є найменшою серед невстановлених вершин. Ця вершина переноситься в масив встановлених вершин і її вихідні ребра релаксуються. Ітерації повторюються, поки всі вершини не будуть перенесені в масив встановлених вершин [2].
Введемо наступні додаткові позначення для формалізації послідовної версії алгоритму Дейкстри: Q – множина нерозглянутих (невстановлених) вершин;
ExtractMin(Q, q,u) – оператор знаходження вершини u такої, що не була розглянута раніше та відстань від початкової вершини до якої ( dqu ) є найменшою: Алгоритм Дейкстри ( D(q) ) можна формалізувати наступним чином:
ExtractMin (Q, q,u) = (mv = d q,Q0 ) × (mi = 0)×i∈Q {(v = Qi )×mv<dqv (mv = d qv ×
×(mi = i) × (+ + i) ∨ E) × (v = Qi ) × erase(Q, mi)}.
SubD1(Q) = ExtractMin(Q, q,u)×v<V {RelaxEdge(q,u,v) × (+ + v)},
D(q) = init(q)×Q≠∅ {SubD1(Q)}.
(
to(D, i)= j<V {setWay (i, j, D j + h j − hi ) × (+ + j)} , де {D j}, j = 0,V – множина найкоротших шляхів з вершини i до всіх інших вершин графу (отриманих алгоритмом Дейкстри).
Використовуючи позначення, можна сформувати послідовну регулярну схему алгоритму Джонсона: α = BellmanFord (V ) ,
SubJohns 1(V )=i<V {to(D(i), i) × (+ + i)} , SubJohns 2(V )=α (UpdateWeig hts × SubJohns1(V ) ∨ E) ,
Johnson = extendGrap h(V ) × SubJohns 2(V ) .
З точки зору оптимізації і прискорення алгоритму, основний акцент має ставитися на прискорення алгоритму Дейкстри, оскільки саме в цьому виконуються основні обчислення. Створення паралельних схем
Класичний підхід до розпаралелювання алгоритмів APSP заключається в паралельному виконанні багатьох ітерацій алгоритму SSSP. Такий підхід є доцільним для виконання обчислень на CPU. При використанні технології GPGPU та моделі CUDA велику роль у підвищенні продуктивності має рівень різноманітності гілок виконання потоків. За моделлю SIMT різні потоки можуть мати різні шляхи виконання інструкцій. Втім, на апаратному рівні різні шляхи виконання виконуються послідовно, навіть якщо певна умова розгалуження виконується лише для одного потоку. Розробники CPU архітектур доклали багато зусиль для створення компенсаційних заходів, зокрема, спекулятивного виконання та передбачення переходів. На GPU таких оптимізацій не було виконано, тому використання розгалужень, шляхи виконання яких сильно відрізняються, призводить до суттєвої втрати продуктивності [7].
В роботі запропоновано паралельні схеми, що ґрунтуються на паралельному виконанні обчислень всередині ітерацій SSSP. Це дає можливість звести рівень розгалуження гілок до мінімуму.
Паралельні схеми алгоритму Дейкстри. Алгоритм Дейкстри складається з двох циклів – зовнішній на кожній ітерації обирає вершину для релаксації всіх ребер, що виходять з неї; внутрішній цикл – власне цикл релаксації ребер. Паралелізація зовнішнього циклу означає, що на кожній ітерації обирається більш ніж одна вершина, ребра якої можуть бути релаксовані без впливу на коректність результатів алгоритму. Паралелізація внутрішнього циклу означає, що вихідні ребра обраної вершини релаксуються одночасно.
Паралелізація зовнішнього циклу. На кожній ітерації алгоритм Дейкстри виконує пошук вершини u такої, що не була розглянута раніше та відстань від початкової вершини до якої ( dqu ) є найменшою. Чим більша кількість обраних вершин на кожній ітерації – тим більший рівень паралелізму.
Таким чином, якщо при послідовній версії обиралася одна вершина u (з умови мінімальності dqu ), то в паралельній версії u – множина вершин v, для яких справджується нерівність: dqv ≤ dqu + Δ , (12) (13) (14) (15)
PExtractMin(Q, q,u) = ExtractLimit(Q, q, m) × size(Q)
• × ∨ (
i=1 dqi<dqm+Δ (add (u,l) ∨ E) × BQ (size(Q)))× eraseAll(Q,u). thrust::device_vector<int> Q(V); INIT<<<blocks, threads>>>(Q, M, C, q); while(Q.size() != 0) { }
RELAX_WITH_MASK<<<blocks, threads>>>(G, M, C); threshold = ExtractMin(Q, C);
UPDATE_MASK<<<blocks, threads>>>(Q, C, M, threshold);
де Δ – довжина найменшого ребра в графі [8]. Такий варіант було обрано для економії пам’яті та
обчислювальної потужності за рахунок малого зниження продуктивності. Для збільшення швидкості виконання алгоритму
(за разунок додаткової пам’яті) можна замінити Δ на Δv , де Δv – довжина найменшого ребра, що виходить з
вершини v. Така заміна дасть можливість корегувати умову «прийнятності» для кожної окремої вершини, що
дозволить збільшити рівень паралелізму. Тому оператор (
ExtractMin (Q, q,u) = ExtractLim it(Q, q, m)×i<size(Q) {dql <dqm +Δ (add (u,l) ∨ E) ×
×(+ + i)}× eraseAll(Q,u).
ExtractLimit – оператор знаходження вершини, відстань від початкової вершини до якої є найменшою
(аналог (
• PSubD1(u) = l∨=1 (v<V (RelaxEdge (q,ul ,v) × (+ + v)} × Bu (size(u))),
Dijkstra(q) = init(q)×Q≠∅ {PExtractMin(Q, q,u) × PSubD1(u)}.
Як видно з формалізованої схеми, кількість потоків дорівнює кількості вершин у графі. Структура алгоритму наведена далі: На кожній ітерації релаксація вихідних ребер відбувається не для всіх вершин, а з певною маскою: 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<int>
thrust::detail::normal_iterator<device_ptr<int>> 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].
Паралелізація внутрішнього циклу. Як було сказано вище, паралелізація внутрішнього циклу означає
одночасну релаксацію всіх ребер обраної вершини. Таким чином, оператор (13) залишається такою ж, як і в
послідовній версії ((
V • PSubD2(u) = ∨ (RelaxEdge(q,u,v) × Bv (v)),
v=1
Dijkstra (q) = init (q)×Q≠∅ {ExtractMin (Q,u) × PSubD 2(u)}.
(16)
Паралелізація алгоритму Беллмана – Форда. Оскільки алгоритм Беллмана – Форда виконується в
алгоритмі Джонсона лише раз, навіть при ідеальній паралельній версії великого виграшу отримати не вдасться.
Втім, алгоритм може опрацьовувати кожну вершину паралельно, тому оператор (
V • RelaxEdges (q) = u∨=1(u<V {RelaxEdge (q,u, v) × (+ + v) × Bv (u)}). (17) Схема всього алгоритму залишається без змін. Аналіз отриманих результатів