Rozwiązywanie problemu najkrótszych ścieżek dla wszystkich par za pomocą algorytmu Floyda-Warshalla w języku C#

Wszyscy rozwiązujemy problem „ ” wiele razy dziennie. Oczywiście nieumyślnie. Rozwiązujemy go w drodze do pracy, przeglądając Internet lub porządkując rzeczy na biurku.

Brzmi trochę… za dużo? Przekonajmy się.

Wyobraź sobie, że postanowiłeś spotkać się ze znajomymi, no cóż... powiedzmy w kawiarni. Przede wszystkim musisz znaleźć trasę (lub ścieżkę) do kawiarni, więc zaczynasz szukać dostępnego transportu publicznego (jeśli idziesz pieszo) lub tras i ulic (jeśli jedziesz samochodem). Szukasz trasy z obecnej lokalizacji do kawiarni, ale nie „żadnej” trasy – najkrótszej lub najszybszej.

Oto jeszcze jeden przykład, który nie jest tak oczywisty jak poprzedni. Podczas pracy nad drogą decydujesz się na „skrót” przez boczną ulicę, ponieważ cóż… to jest „skrót” i „szybciej” jest iść tą drogą. Ale skąd wiesz, że ten „skrót” jest szybszy? Na podstawie osobistych doświadczeń rozwiązałeś problem „najkrótszej ścieżki” i wybrałeś trasę, która prowadzi przez boczną ulicę.

W obu przykładach „najkrótsza” trasa jest określana albo odległością, albo czasem potrzebnym na dotarcie z jednego miejsca do drugiego. Przykłady podróży są bardzo naturalnymi zastosowaniami i w szczególności problemu „najkrótszej ścieżki”. Co jednak z przykładem układania rzeczy na biurku? W tym przypadku „najkrótsza” może reprezentować liczbę lub złożoność działań, które musisz wykonać, aby na przykład uzyskać kartkę papieru: otwórz biurko, otwórz folder, weź kartkę papieru vs. weź kartkę papieru prosto z biurka.

Wszystkie powyższe przykłady przedstawiają problem znalezienia najkrótszej ścieżki między dwoma wierzchołkami w grafie (trasa lub ścieżka między dwoma miejscami, liczba działań lub złożoność pobrania kartki papieru z jednego miejsca lub drugiego). Ta klasa problemów najkrótszej ścieżki nazywana jest SSSP (Single Source Shortest Path) , a podstawowym algorytmem rozwiązywania tych problemów jest , który ma złożoność obliczeniową O(n^2) .

Ale czasami musimy znaleźć wszystkie najkrótsze ścieżki pomiędzy wszystkimi wierzchołkami. Rozważmy następujący przykład: tworzysz mapę dla swoich regularnych ruchów pomiędzy domem , pracą i teatrem . W tym scenariuszu otrzymasz 6 tras: work ⭢ home , home ⭢ work , work ⭢ theatre , theatre ⭢ work , home ⭢ theatre i theatre ⭢ home (trasy powrotne mogą być różne ze względu na drogi jednokierunkowe na przykład).

Dodanie większej ilości miejsc na mapie spowoduje znaczący wzrost liczby kombinacji – zgodnie z można to obliczyć następująco:

 A(k, n) = n! / (n - m)! // where n - is a number of elements, // k - is a number of elements in permutation (in our case k = 2)

Co daje nam 12 tras dla 4 miejsc i 90 tras dla 10 miejsc (co jest imponujące). Uwaga… to nie uwzględnia punktów pośrednich między miejscami, np. aby dostać się z domu do pracy musisz przejść przez 4 ulice, iść wzdłuż rzeki i przejść przez most. Jeśli sobie wyobrazimy, niektóre trasy mogą mieć wspólne punkty pośrednie… cóż… w rezultacie będziemy mieli bardzo duży graf z wieloma wierzchołkami, gdzie każdy wierzchołek będzie reprezentował albo miejsce, albo znaczący punkt pośredni. Klasa problemów, w których musimy znaleźć wszystkie najkrótsze ścieżki między wszystkimi parami wierzchołków w grafie, nazywa się APSP (All Pairs Shortest Paths), a podstawowym algorytmem rozwiązywania tych problemów jest , który ma złożoność obliczeniową O(n^3) .

I to jest algorytm, który dziś wdrożymy 🙂

Algorytm Floyda-Warshalla

Algorytm Floyda-Warshalla znajduje wszystkie najkrótsze ścieżki między każdą parą wierzchołków w grafie. Algorytm został opublikowany przez Roberta Floyda w [1] (więcej szczegółów można znaleźć w sekcji „Odniesienia”). W tym samym roku Peter Ingerman w [2] opisał nowoczesną implementację algorytmu w formie trzech zagnieżdżonych pętli for :

 algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm // where W - is a weight matrix of N x N size, // min() - is a function which returns lesser of it's arguments

Jeśli nigdy nie miałeś okazji pracować z grafem przedstawionym w formie macierzy, może być trudno zrozumieć, co robi powyższy algorytm. Więc, aby upewnić się, że jesteśmy na tej samej stronie, przyjrzyjmy się, jak graf może być przedstawiony w formie macierzy i dlaczego taka reprezentacja jest korzystna dla rozwiązania problemu najkrótszej ścieżki.

Poniższy rysunek ilustruje graf 5 wierzchołków. Po lewej stronie graf jest przedstawiony w formie wizualnej, która składa się z okręgów i strzałek, gdzie każdy okrąg reprezentuje wierzchołek, a strzałka reprezentuje krawędź z kierunkiem. Liczba wewnątrz okręgu odpowiada numerowi wierzchołka, a liczba nad krawędzią odpowiada wadze krawędzi. Po prawej stronie ten sam graf jest przedstawiony w formie macierzy wag. Macierz wag jest formą , w której każda komórka macierzy zawiera „wagę” – odległość między wierzchołkiem i (wiersz) a wierzchołkiem j (kolumna). Macierz wag nie zawiera informacji o „ścieżce” między wierzchołkami (lista wierzchołków, przez którą przechodzisz z i do j ) – tylko wagę (odległość) między tymi wierzchołkami.

W macierzy wag możemy zobaczyć, że wartości komórek są równe wagom między wierzchołkami. Dlatego też, jeśli zbadamy ścieżki od wierzchołka 0 (wiersz 0 ), zobaczymy, że… istnieje tylko jedna ścieżka – od 0 do 1 Jednak na reprezentacji wizualnej możemy wyraźnie zobaczyć ścieżki od wierzchołka 0 do wierzchołków 2 i 3 (przez wierzchołek 1 ). Powód tego jest prosty – w stanie początkowym macierz wag zawiera odległość tylko między sąsiadującymi wierzchołkami. Jednak sama ta informacja wystarcza, aby znaleźć resztę.

Zobaczmy jak to działa. Zwróć uwagę na komórkę W[0, 1] . Jej wartość wskazuje, że istnieje ścieżka od wierzchołka 0 do wierzchołka 1 o wadze równej 1 (w skrócie możemy zapisać to jako: 0 ⭢ 1: 1 ). Mając tę wiedzę, możemy teraz przeskanować wszystkie komórki wiersza 1 (który zawiera wszystkie wagi wszystkich ścieżek z wierzchołka 1 ) i przenieść te informacje z powrotem do wiersza 0 , zwiększając wagę o 1 (wartość W[0, 1] ).

Używając tych samych kroków możemy znaleźć ścieżki od wierzchołka 0 przez inne wierzchołki. Podczas wyszukiwania może się zdarzyć, że istnieje więcej niż jedna ścieżka prowadząca do tego samego wierzchołka i co ważniejsze wagi tych ścieżek mogą być różne. Przykład takiej sytuacji jest zilustrowany na poniższym rysunku, gdzie wyszukiwanie od wierzchołka 0 przez wierzchołek 2 ujawniło jeszcze jedną ścieżkę do wierzchołka 3 o mniejszej wadze.

Mamy dwie ścieżki: oryginalną ścieżkę – 0 ⭢ 3: 4 i nową ścieżkę, którą właśnie odkryliśmy – 0 ⭢ 2 ⭢ 3: 3 (pamiętaj, macierz wag nie zawiera ścieżek, więc nie wiemy, które wierzchołki są zawarte w oryginalnej ścieżce). To moment, w którym podejmujemy decyzję o zachowaniu najkrótszej ścieżki i zapisaniu 3 w komórce W[0, 3] .

Wygląda na to, że właśnie znaleźliśmy naszą pierwszą najkrótszą ścieżkę!

Kroki, które właśnie widzieliśmy, są istotą algorytmu Floyda-Warshalla. Spójrzmy jeszcze raz na pseudokod algorytmu:

 algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm

Najbardziej zewnętrzny cykl for on k iteruje po wszystkich wierzchołkach w grafie i w każdej iteracji zmienna k reprezentuje wierzchołek , przez który przeszukujemy ścieżki . Wewnętrzny cykl for on i również iteruje po wszystkich wierzchołkach w grafie i w każdej iteracji i reprezentuje wierzchołek , przez który przeszukujemy ścieżki . I na koniec najbardziej wewnętrzny cykl for on j iteruje po wszystkich wierzchołkach w grafie i w każdej iteracji j reprezentuje wierzchołek , przez który przeszukujemy ścieżki. W połączeniu daje nam to następujące: w każdej iteracji k przeszukujemy ścieżki ze wszystkich wierzchołków i do wszystkich wierzchołków j przez wierzchołek k . Wewnątrz cyklu porównujemy ścieżkę i ⭢ j (reprezentowaną przez W[i, j] ) ze ścieżką i ⭢ k ⭢ j (reprezentowaną przez sumę W[I, k] i W[k, j] ) i zapisujemy najkrótszą z powrotem do W[i, j] .

Teraz, gdy rozumiemy już mechanikę, czas na wdrożenie algorytmu.

Podstawowa implementacja

Zanim przejdziemy do realizacji, musimy wyjaśnić kilka kwestii technicznych:

1. Wszystkie implementacje działają z macierzą wag reprezentowaną w formie tablicy liniowej.
2. Wszystkie implementacje używają arytmetyki liczb całkowitych. Brak ścieżki między wierzchołkami jest reprezentowany przez specjalną stałą wagę: NO_EDGE = (int.MaxValue / 2) - 1 .

Teraz zastanówmy się, dlaczego tak jest.

Odnośnie #1. Kiedy mówimy o macierzach, definiujemy komórki w kategoriach „wierszy” i „kolumn”. Z tego powodu wydaje się naturalne wyobrażanie sobie macierzy w formie „kwadratu” lub „prostokąta” (Rysunek 4a).

Jednakże nie oznacza to koniecznie, że musimy przedstawić macierz w formie tablicy tablic (Rysunek 4b), aby trzymać się naszej wyobraźni. Zamiast tego możemy przedstawić macierz w formie tablicy liniowej (Rysunek 4c), gdzie indeks każdej komórki jest obliczany przy użyciu następującego wzoru:

 i = row * row_size + col; // where row - cell row index, // col - cell column index, // row_size - number of cells in a row.

Liniowa tablica macierzy wag jest warunkiem koniecznym skutecznego wykonania algorytmu Floyda-Warshalla. Powód tego nie jest prosty i szczegółowe wyjaśnienie zasługuje na osobny post… lub kilka postów. Jednak obecnie ważne jest, aby wspomnieć, że taka reprezentacja znacząco poprawia , co w rzeczywistości ma duży wpływ na wydajność algorytmu.

Odnośnie #2. Jeśli przyjrzysz się bliżej pseudokodowi algorytmu, nie znajdziesz żadnych sprawdzeń związanych z istnieniem ścieżki między dwoma wierzchołkami. Zamiast tego pseudokod po prostu używa funkcji min() . Powód jest prosty – pierwotnie, jeśli nie ma ścieżki między dwoma wierzchołkami, wartość komórki jest ustawiana na infinity i we wszystkich językach, z wyjątkiem JavaScript, wszystkie wartości są mniejsze od infinity . W przypadku liczb całkowitych może być kuszące użycie int.MaxValue jako wartości „brak ścieżki”. Jednak doprowadzi to do przepełnienia całkowitego w przypadkach, gdy wartości obu ścieżek i ⭢ k i k ⭢ j będą równe int.MaxValue . Dlatego używamy wartości, która jest o jeden mniejsza niż połowa int.MaxValue .

Jest to rzeczywiście możliwe, ale niestety doprowadzi to do znacznego spadku wydajności. Krótko mówiąc, CPU przechowuje statystyki wyników oceny rozgałęzień, np. gdy niektóre z instrukcji if są oceniane jako true lub false . Wykorzystuje te statystyki do wykonania kodu „statystycznie przewidzianej rozgałęzienia” z góry, zanim zostanie oceniona rzeczywista instrukcja if (nazywa się to ), a zatem wykonuje kod wydajniej. Jednak gdy przewidywanie CPU jest niedokładne, powoduje to znaczną utratę wydajności w porównaniu z prawidłowym przewidywaniem i wykonaniem bezwarunkowym, ponieważ CPU musi zatrzymać się i obliczyć prawidłową rozgałęzienie.

Ponieważ w każdej iteracji k aktualizujemy znaczną część macierzy wag, statystyki rozgałęzień procesora stają się bezużyteczne, ponieważ nie ma wzorca kodu, wszystkie rozgałęzienia są oparte wyłącznie na danych. Tak więc takie sprawdzenie spowoduje znaczną liczbę błędnych .

Tutaj również proszę Was, abyście mi uwierzyli (na razie), a potem poświęcili trochę czasu na zgłębienie tematu.

Uhh, wygląda na to, że skończyliśmy część teoretyczną – zaimplementujmy algorytm (oznaczmy tę implementację jako Baseline ):

 public void Baseline(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

Powyższy kod jest niemal identyczną kopią wcześniej wspomnianego pseudokodu, z jednym wyjątkiem – zamiast Math.Min() używamy if aby zaktualizować komórkę tylko wtedy, gdy jest to konieczne.

Powód jest prosty. W chwili pisania tego tekstu JIT emituje prawie równoważny kod dla obu implementacji if i Math.Min . Możesz to sprawdzić szczegółowo na , ale oto fragmenty ciał pętli głównej:

 // // == Assembly code of implementation of innermost loop for of j using if. // 53: movsxd r14, r14d // // compare matrix[i * sz + j] and distance (Condition) // 56: cmp [rdx+r14*4+0x10], ebx 5b: jle short 64 // // if matrix[i * sz + j] greater than distance write distance to matrix // 5d: movsxd rbp, ebp 60: mov [rdx+rbp*4+0x10], ebx 64: // end of loop // // == Assembly code of implementation of innermost loop for of j using Math.Min. // 4f: movsxd rbp, ebp 52: mov r14d, [rdx+rbp*4+0x10] // // compare matrix[i * sz + j] and distance (Condition) // 57: cmp r14d, ebx. // // if matrix[i * sz + j] less than distance write value to matrix // 5a: jle short 5e // // otherwise write distance to matrix // 5c: jmp short 61 5e: mov ebx, r14d 61: mov [rdx+rbp*4+0x10], ebx 65: // end of loop

Możesz zobaczyć, że niezależnie od tego, czy używamy if czy Math.Min nadal istnieje warunkowe sprawdzenie. Jednak w przypadku if nie ma zbędnego zapisu.

Teraz, gdy zakończyliśmy implementację, czas uruchomić kod i zobaczyć, jak szybki jest?! Używam (wersja 0.12.1) do testowania kodu. Wszystkie grafy używane w testach porównawczych są grafami o 300, 600, 1200, 2400 i 4800 wierzchołkach. Liczba krawędzi w grafach wynosi około 80% możliwego maksimum, co dla skierowanych, acyklicznych grafów można obliczyć jako:

 var max = v * (v - 1)) / 2; // where v - is a number of vertexes in a graph.

Dajmy czadu!

Oto wyniki testów porównawczych przeprowadzonych na moim komputerze (Windows 10.0.19042, procesor Intel Core i7-7700 3,60 GHz (Kaby Lake) / 8 procesorów logicznych / 4 rdzenie):

Metoda	Rozmiar	Mieć na myśli	Błąd	Odchylenie standardowe
Linia bazowa	300	27,525 milisekundy	0,1937 milisekundy	0,1617 milisekundy
Linia bazowa	600	217,897 milisekund	1,6415 milisekundy	1,5355 milisekundy
Linia bazowa	1200	1,763.335 milisekund	7,4561 milisekundy	6,2262 milisekundy
Linia bazowa	2400	14 533,335 milisekund	63,3518 milisekund	52,9016 milisekund
Linia bazowa	4800	119 768,219 milisekund	181,5514 milisekundy	160,9406 milisekund

Z wyników możemy wywnioskować, że czas obliczeń rośnie drastycznie w porównaniu do rozmiaru grafu – dla grafu o 300 wierzchołkach zajęło to 27 milisekund, dla grafu o 2400 wierzchołkach – 14,5 sekundy, a dla grafu o 4800 wierzchołkach – 119 sekund, czyli prawie 2 minuty !

Patrząc na kod algorytmu, trudno sobie wyobrazić, że jest coś, co możemy zrobić, aby przyspieszyć obliczenia… ponieważ… są TRZY pętle, tylko TRZY pętle.

Jednak, jak to zwykle bywa – możliwości kryją się w szczegółach 🙂

Poznaj swoje dane – rzadkie wykresy

Algorytm Floyda-Warshalla jest podstawowym algorytmem rozwiązywania problemu najkrótszej ścieżki dla wszystkich par, zwłaszcza w przypadku grafów i (ponieważ algorytm przeszukuje ścieżki pomiędzy wszystkimi parami wierzchołków).

Jednak w naszych eksperymentach używamy skierowanych, acyklicznych grafów, które mają wspaniałą właściwość – jeśli istnieje ścieżka od wierzchołka 1 do wierzchołka 2 , to nie ma ścieżki od wierzchołka 2 do wierzchołka 1 Dla nas oznacza to, że istnieje wiele nieprzylegających wierzchołków, które możemy pominąć, jeśli nie ma ścieżki od i do k (oznaczamy tę implementację jako SpartialOptimisation ).

 public void SpartialOptimisation(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

Poniżej przedstawiono wyniki poprzednich ( Baseline ) i obecnych ( SpartialOptimisation ) implementacji na tym samym zestawie wykresów (najszybsze wyniki wyróżniono pogrubioną czcionką ):

Metoda	Rozmiar	Mieć na myśli	Błąd	Odchylenie standardowe	Stosunek
Linia bazowa	300	27,525 milisekundy	0,1937 milisekundy	0,1617 milisekundy	1,00
Częściowa optymalizacja	300	12,399 milisekund	0,0943 milisekundy	0,0882 milisekundy	0,45
Linia bazowa	600	217,897 milisekund	1,6415 milisekundy	1,5355 milisekundy	1,00
Częściowa optymalizacja	600	99,122 milisekundy	0,8230 milisekundy	0,7698 milisekundy	0,45
Linia bazowa	1200	1,763.335 milisekund	7,4561 milisekundy	6,2262 milisekundy	1,00
Częściowa optymalizacja	1200	766,675 milisekundy	6,1147 milisekund	5,7197 milisekundy	0,43
Linia bazowa	2400	14 533,335 milisekund	63,3518 milisekund	52,9016 milisekund	1,00
Częściowa optymalizacja	2400	6,507.878 milisekund	28,2317 milisekund	26,4079 milisekundy	0,45
Linia bazowa	4800	119 768,219 milisekund	181,5514 milisekundy	160,9406 milisekund	1,00
Częściowa optymalizacja	4800	55 590,374 milisekundy	414,6051 milisekundy	387,8218 milisekund	0,46

Naprawdę imponujące!

Skróciliśmy czas obliczeń o połowę! Oczywiście, im gęstszy wykres, tym mniejsze będzie przyspieszenie. Jednak jest to jedna z optymalizacji, która przydaje się, jeśli wiesz z góry, z jaką klasą danych zamierzasz pracować.

Czy możemy zrobić więcej? 🙂

Poznaj swój sprzęt – paralelizm

---

Mój komputer jest wyposażony w procesor Inter Core i7-7700 CPU 3.60GHz (Kaby Lake) z 8 logicznymi procesorami ( ) lub 4 rdzeniami z technologią . Posiadanie więcej niż jednego rdzenia jest jak posiadanie większej ilości „zapasowych rąk”, które możemy wykorzystać. Zobaczmy więc, która część pracy może być wydajnie podzielona między wielu pracowników, a następnie sparalelizowana.

Pętle są zawsze najbardziej oczywistym kandydatem do paralelizacji, ponieważ w większości przypadków iteracje są niezależne i mogą być wykonywane jednocześnie. W algorytmie mamy trzy pętle, które powinniśmy przeanalizować i sprawdzić, czy istnieją jakieś zależności, które uniemożliwiają nam ich paralelizację.

Zacznijmy od pętli for of k . W każdej iteracji pętla oblicza ścieżki z każdego wierzchołka do każdego wierzchołka przez wierzchołek k . Nowe i zaktualizowane ścieżki są przechowywane w macierzy wag. Patrząc na iteracje możesz zauważyć – mogą być wykonywane w dowolnej kolejności: 0, 1, 2, 3 lub 3, 1, 2, 0 bez wpływu na wynik. Jednak nadal muszą być wykonywane po kolei, w przeciwnym razie niektóre iteracje nie będą mogły używać nowych lub zaktualizowanych ścieżek, ponieważ nie zostaną jeszcze zapisane w macierzy wag. Taka niespójność z pewnością zmiażdży wyniki. Więc musimy dalej szukać.

Następnym kandydatem jest for of i . W każdej iteracji pętla odczytuje ścieżkę z wierzchołka i do wierzchołka k (komórka: W[i, k] ), ścieżkę z wierzchołka k do wierzchołka j (komórka: W[k, j ]), a następnie sprawdza, czy znana ścieżka z i do j (komórka: W[i, j] ) jest krótsza niż ścieżka i ⭢ k ⭢ j (suma: W[i, k] + W[k, j] ). Jeśli przyjrzysz się bliżej wzorcowi dostępu, możesz zauważyć – w każdej iteracji pętla i odczytuje dane z k wiersza i aktualizowanego i wiersza, co zasadniczo oznacza – iteracje są niezależne i mogą być wykonywane nie tylko w dowolnej kolejności, ale także równolegle!

Wygląda to obiecująco, więc zaimplementujmy to (oznaczymy tę implementację jako SpartialParallelOptimisations ).

 public void SpartialParallelOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }

Oto wyniki implementacji Baseline , SpartialOptimisation i SpartialParallelOptimisations na tym samym zestawie wykresów (paralelizacja jest realizowana przy użyciu klasy ):

Metoda	Rozmiar	Mieć na myśli	Błąd	Odchylenie standardowe	Stosunek
Linia bazowa	300	27,525 milisekundy	0,1937 milisekundy	0,1617 milisekundy	1,00
Częściowa optymalizacja	300	12,399 milisekundy	0,0943 milisekundy	0,0882 milisekundy	0,45
Częściowe Równoległe Optymalizacje	300	6,252 milisekundy	0,0211 milisekundy	0,0187 milisekundy	0,23
Linia bazowa	600	217,897 milisekund	1,6415 milisekundy	1,5355 milisekundy	1,00
Częściowa optymalizacja	600	99,122 milisekundy	0,8230 milisekundy	0,7698 milisekundy	0,45
Częściowe Równoległe Optymalizacje	600	35,825 milisekundy	0,1003 milisekundy	0,0837 milisekundy	0,16
Linia bazowa	1200	1,763.335 milisekund	7,4561 milisekundy	6,2262 milisekundy	1,00
Częściowa optymalizacja	1200	766,675 milisekundy	6,1147 milisekund	5,7197 milisekundy	0,43
Częściowe Równoległe Optymalizacje	1200	248,801 milisekundy	0,6040 milisekundy	0,5043 milisekundy	0,14
Linia bazowa	2400	14 533,335 milisekund	63,3518 milisekund	52,9016 milisekund	1,00
Częściowa optymalizacja	2400	6,507.878 milisekund	28,2317 milisekund	26,4079 milisekundy	0,45
Częściowe Równoległe Optymalizacje	2400	2,076.403 milisekundy	20,8320 milisekund	19,4863 milisekundy	0,14
Linia bazowa	4800	119 768,219 milisekund	181,5514 milisekundy	160,9406 milisekund	1,00
Częściowa optymalizacja	4800	55 590,374 milisekundy	414,6051 milisekundy	387,8218 milisekund	0,46
Częściowe Równoległe Optymalizacje	4800	15 614,506 milisekund	115,6996 milisekund	102,5647 milisekundy	0,13

Z wyników możemy wywnioskować, że paralelizacja pętli for of i skróciła czas obliczeń o 2-4 razy w porównaniu do poprzedniej implementacji ( SpartialOptimisation )! To bardzo imponujące, jednak należy pamiętać, że paralelizacja czystych obliczeń zużywa wszystkie dostępne zasoby obliczeniowe, co może prowadzić do niedoboru zasobów innych aplikacji w systemie. Paralelizację należy stosować ostrożnie.

Jak pewnie zgadliście – to nie koniec 🙂

Poznaj swoją platformę – wektoryzacja

to przekształcenie kodu operującego na pojedynczym elemencie w kod operujący na wielu elementach jednocześnie. Może się to wydawać skomplikowane, ale zobaczmy, jak to działa na prostym przykładzie:

 var a = new [] { 5, 7, 8, 1 }; var b = new [] { 4, 2, 2, 6 }; var c = new [] { 0, 0, 0, 0 }; for (var i = 0; i < 4; ++i) c[i] = a[i] + b[i];

W uproszczeniu wykonanie iteracji 0 powyższej pętli for na poziomie procesora można zilustrować następująco:

Oto, co się dzieje. CPU ładuje wartości tablic a i b z pamięci do rejestrów CPU (jeden rejestr może przechowywać dokładnie jedną wartość). Następnie CPU wykonuje operację dodawania skalarnego na tych rejestrach i zapisuje wynik z powrotem do pamięci głównej – bezpośrednio do tablicy c . Ten proces jest powtarzany cztery razy, zanim pętla się zakończy.

Wektoryzacja oznacza wykorzystanie specjalnych rejestrów procesora – rejestrów wektorowych lub SIMD (single instruction multiple data) i odpowiadających im instrukcji procesora w celu wykonywania operacji na wielu wartościach tablicy jednocześnie:

Oto, co się dzieje. CPU ładuje wartości tablic a i b z pamięci do rejestrów CPU (jednak tym razem jeden rejestr wektorowy może przechowywać dwie wartości tablicowe). Następnie CPU wykonuje operację dodawania wektorów na tych rejestrach i zapisuje wynik z powrotem do pamięci głównej – bezpośrednio do tablicy c . Ponieważ operujemy na dwóch wartościach jednocześnie, proces ten powtarza się dwa razy zamiast czterech.

Aby zaimplementować wektoryzację w .NET możemy użyć typów i (możemy również użyć typów z przestrzeni nazw , jednak są one nieco zaawansowane dla obecnej implementacji, więc nie będę ich używał, ale możesz spróbować samodzielnie):

 public void SpartialVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } }

Wektoryzowany kod może wyglądać nieco dziwnie, więc przeanalizujmy go krok po kroku.

Tworzymy wektor ścieżki i ⭢ k : var ik_vec = new Vector<int>(matrix[i * sz + k]) . W rezultacie, jeśli wektor może pomieścić cztery wartości typu int , a waga ścieżki i ⭢ k jest równa 5, utworzymy wektor czterech 5 – [5, 5, 5, 5]. Następnie, w każdej iteracji, jednocześnie obliczamy ścieżki od wierzchołka i do wierzchołków j, j + 1, ..., j + Vector<int>.Count .

Cały proces obliczeniowy można podzielić na trzy etapy:

1. Załaduj informacje z macierzy wag do wektorów: ij_vec i ikj_vec .
2. Porównaj wektory ij_vec i ikj_vec , a następnie wybierz najmniejsze wartości do nowego wektora r_vec .
3. Zapisz r_vec z powrotem do macierzy wag.

Podczas gdy #1 i #3 są dość proste, #2 wymaga wyjaśnienia. Jak wspomniano wcześniej, w przypadku wektorów manipulujemy wieloma wartościami jednocześnie. Dlatego wynik niektórych operacji nie może być osobliwy, np. operacja porównania Vector.LessThan(ij_vec, ikj_vec) nie może zwrócić true lub false , ponieważ porównuje wiele wartości. Zamiast tego zwraca nowy wektor, który zawiera wyniki porównania między odpowiadającymi sobie wartościami z wektorów ij_vec i ikj_vec ( -1 , jeśli wartość z ij_vec jest mniejsza niż wartość z ikj_vec i 0 jeśli jest inaczej). Zwrócony wektor (sam w sobie) nie jest zbyt przydatny, jednak możemy użyć operacji wektorowej Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec) , aby wyodrębnić wymagane wartości z wektorów ij_vec i ikj_vec do nowego wektora – r_vec . Operacja ta zwraca nowy wektor, w którym wartości są równe mniejszej z dwóch odpowiadających sobie wartości wektorów wejściowych, tj. jeśli wartość wektora lt_vec jest równa -1 , to operacja wybiera wartość z ij_vec , w przeciwnym wypadku wybiera wartość z ikj_vec .

Oprócz #2 , jest jeszcze jedna część, która wymaga wyjaśnienia – druga, niewektoryzowana pętla. Jest ona wymagana, gdy rozmiar macierzy wag nie jest iloczynem wartości Vector<int>.Count . W takim przypadku główna pętla nie może przetworzyć wszystkich wartości (ponieważ nie można załadować części wektora), a druga, niewektoryzowana pętla obliczy pozostałe iteracje.

Oto rezultaty tej i wszystkich poprzednich implementacji:

Metoda	sz	Mieć na myśli	Błąd	Odchylenie standardowe	Stosunek
Linia bazowa	300	27,525 milisekundy	0,1937 milisekundy	0,1617 milisekundy	1,00
Częściowa optymalizacja	300	12,399 milisekund	0,0943 milisekundy	0,0882 milisekundy	0,45
Częściowe Równoległe Optymalizacje	300	6,252 milisekundy	0,0211 milisekundy	0,0187 milisekundy	0,23
SpartialVectorOptimisations	300	3,056 milisekundy	0,0301 milisekundy	0,0281 milisekundy	0,11
Linia bazowa	600	217,897 milisekund	1,6415 milisekundy	1,5355 milisekundy	1,00
Częściowa optymalizacja	600	99,122 milisekundy	0,8230 milisekundy	0,7698 milisekundy	0,45
Częściowe Równoległe Optymalizacje	600	35,825 milisekundy	0,1003 milisekundy	0,0837 milisekundy	0,16
SpartialVectorOptimisations	600	24,378 milisekundy	0,4320 milisekundy	0,4041 milisekundy	0,11
Linia bazowa	1200	1,763.335 milisekund	7,4561 milisekundy	6,2262 milisekundy	1,00
Częściowa optymalizacja	1200	766,675 milisekundy	6,1147 milisekund	5,7197 milisekundy	0,43
Częściowe Równoległe Optymalizacje	1200	248,801 milisekundy	0,6040 milisekundy	0,5043 milisekundy	0,14
SpartialVectorOptimisations	1200	185,628 milisekund	2,1240 milisekundy	1,9868 milisekundy	0,11
Linia bazowa	2400	14 533,335 milisekund	63,3518 milisekund	52,9016 milisekund	1,00
Częściowa optymalizacja	2400	6,507.878 milisekund	28,2317 milisekund	26,4079 milisekundy	0,45
Częściowe Równoległe Optymalizacje	2400	2,076.403 milisekundy	20,8320 milisekund	19,4863 milisekundy	0,14
SpartialVectorOptimisations	2400	2,568,676 milisekund	31,7359 milisekund	29,6858 milisekund	0,18
Linia bazowa	4800	119 768,219 milisekund	181,5514 milisekund	160,9406 milisekund	1,00
Częściowa optymalizacja	4800	55 590,374 milisekundy	414,6051 milisekundy	387,8218 milisekund	0,46
Częściowe Równoległe Optymalizacje	4800	15 614,506 milisekund	115,6996 milisekund	102,5647 milisekundy	0,13
SpartialVectorOptimisations	4800	18 257,991 milisekund	84,5978 milisekund	79,1329 milisekundy	0,15

Z wyniku wynika, że wektoryzacja znacznie skróciła czas obliczeń – od 3 do 4 razy w porównaniu do wersji nieparalelizowanej ( SpartialOptimisation ). Interesującym momentem jest to, że wersja zwektoryzowana przewyższa wersję równoległą ( SpartialParallelOptimisations ) również na mniejszych grafach (mniej niż 2400 wierzchołków).

I na koniec – połączmy wektoryzację i paralelizm!

Poznaj swoją platformę i sprzęt – wektoryzację i paralelizm!

Ostatnia implementacja łączy w sobie wysiłki paralelizacji i wektoryzacji i… szczerze mówiąc, brakuje jej indywidualności 🙂 Ponieważ… cóż, właśnie zastąpiliśmy ciało Parallel.For ze SpartialParallelOptimisations zwektoryzowaną pętlą ze SpartialVectorOptimisations :

 public void SpartialParallelVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } }

Wszystkie wyniki tego posta są przedstawione poniżej. Zgodnie z oczekiwaniami, jednoczesne użycie paralelizmu i wektoryzacji wykazało najlepsze wyniki, skracając czas obliczeń nawet 25 razy (dla grafów 1200 wierzchołków) w porównaniu z początkową implementacją.

Metoda	sz	Mieć na myśli	Błąd	Odchylenie standardowe	Stosunek
Linia bazowa	300	27,525 milisekundy	0,1937 milisekundy	0,1617 milisekundy	1,00
Częściowa optymalizacja	300	12,399 milisekundy	0,0943 milisekundy	0,0882 milisekundy	0,45
Częściowe Równoległe Optymalizacje	300	6,252 milisekundy	0,0211 milisekundy	0,0187 milisekundy	0,23
SpartialVectorOptimisations	300	3,056 milisekundy	0,0301 milisekundy	0,0281 milisekundy	0,11
Częściowe Optymalizacje Wektorów Równoległych	300	3,008 milisekundy	0,0075 milisekundy	0,0066 milisekundy	0,11
Linia bazowa	600	217,897 milisekund	1,6415 milisekundy	1,5355 milisekundy	1,00
Częściowa optymalizacja	600	99,122 milisekundy	0,8230 milisekundy	0,7698 milisekundy	0,45
Częściowe Równoległe Optymalizacje	600	35,825 milisekundy	0,1003 milisekundy	0,0837 milisekundy	0,16
SpartialVectorOptimisations	600	24,378 milisekundy	0,4320 milisekundy	0,4041 milisekundy	0,11
Częściowe Optymalizacje Wektorów Równoległych	600	13,425 milisekundy	0,0319 milisekundy	0,0283 milisekundy	0,06
Linia bazowa	1200	1,763.335 milisekund	7,4561 milisekundy	6,2262 milisekundy	1,00
Częściowa optymalizacja	1200	766,675 milisekundy	6,1147 milisekund	5,7197 milisekundy	0,43
Częściowe Równoległe Optymalizacje	1200	248,801 milisekundy	0,6040 milisekundy	0,5043 milisekundy	0,14
SpartialVectorOptimisations	1200	185,628 milisekund	2,1240 milisekundy	1,9868 milisekundy	0,11
Częściowe Optymalizacje Wektorów Równoległych	1200	76,770 milisekund	0,3021 milisekundy	0,2522 milisekundy	0,04
Linia bazowa	2400	14 533,335 milisekund	63,3518 milisekund	52,9016 milisekund	1,00
Częściowa optymalizacja	2400	6,507.878 milisekund	28,2317 milisekund	26,4079 milisekundy	0,45
Częściowe Równoległe Optymalizacje	2400	2,076.403 milisekundy	20,8320 milisekund	19,4863 milisekundy	0,14
SpartialVectorOptimisations	2400	2,568,676 milisekund	31,7359 milisekund	29,6858 milisekund	0,18
Częściowe Optymalizacje Wektorów Równoległych	2400	1,281.877 milisekund	25,1127 milisekundy	64,8239 milisekund	0,09
Linia bazowa	4800	119 768,219 milisekund	181,5514 milisekundy	160,9406 milisekund	1,00
Częściowa optymalizacja	4800	55 590,374 milisekundy	414,6051 milisekundy	387,8218 milisekund	0,46
Częściowe Równoległe Optymalizacje	4800	15 614,506 milisekund	115,6996 milisekund	102,5647 milisekundy	0,13
SpartialVectorOptimisations	4800	18 257,991 milisekund	84,5978 milisekund	79,1329 milisekundy	0,15
Częściowe Optymalizacje Wektorów Równoległych	4800	12 785,824 milisekund	98,6947 milisekund	87,4903 milisekundy	0,11

Wniosek

W tym poście nieco zagłębiliśmy się w problem najkrótszej ścieżki dla wszystkich par i zaimplementowaliśmy algorytm Floyda-Warshalla w C#, aby go rozwiązać. Zaktualizowaliśmy również naszą implementację, aby honorować dane i wykorzystywać funkcje niskiego poziomu .NET i sprzętu.

W tym poście zagraliśmy „all in”. Jednak w rzeczywistych zastosowaniach ważne jest, aby pamiętać – nie jesteśmy sami. Hardkorowa paralelizacja może znacznie obniżyć wydajność systemu i spowodować więcej szkody niż pożytku. Z drugiej strony wektoryzacja jest nieco bezpieczniejsza, jeśli jest wykonywana w sposób niezależny od platformy. Pamiętaj, że niektóre instrukcje wektorowe mogą być dostępne tylko na niektórych procesorach.

Mam nadzieję, że lektura przypadła Ci do gustu i że dobrze się bawiłeś, „wyciskając” trochę wydajności z algorytmu zaledwie w TRZECH cyklach 🙂

Odniesienia

---