Zaawansowane algorytmy i struktury danych/Wykład 13

From Studia Informatyczne

Algorytmy równoległe I

Spis treści


W module tym zajmiemy się przyspieszaniem obliczeń za pomocą korzystania z wielu procesorów (maszyn) działających równolegle. Niestety nie ma ogólnie przyjętego modelu obliczeń równoległych. Rozważymy w tym module dwa modele: maszynę PRAM i układy arytmetyczne (logiczne). O ile maszyna PRAM jest modelem wysoko-poziomowym, to układy arytmetyczne są modelem niskopoziomowym, ale niewątpliwie bardzo istotnym.

Model równoległej abstrakcyjnej Maszyny PRAM

Na początku rozważymy wyidealizowany model obliczeń równoległych zwany Równoległą Maszyną ze Swobodnym Dostępem do Pamięci, w skrócie PRAM (od ang. Parallel Random Access Machine, wymawiany piram).

Grafika:Parallel1-1.png
Rysunek 1: Struktura koncepcyjna PRAMu

Maszyna PRAM składa się z wielu procesorów pracujących synchronicznie, korzystących ze wspólnej pamięci (która oprócz przechowywania danych służy do komunikacji między procesorami). Każdy procesor jest standardowym komputerem typu RAM (ang. Random Access Machine). Zakładamy, że procesory są ponumerowane liczbami naturalnymi. Procesory wykonują jeden wspólny program, ale wykonanie poszczególnych instrukcji zależy od indeksu procesora. W jednym kroku procesor pobiera dane z pamięci, potem wykonuje operację, którą może być wpisanie pewnych danych. Wszystkie procesory wykonują jeden krok jednocześnie. Równoległość jest wyrażona poprzez następującą instrukcję:

\textbf{for all } i \in X  \textbf{do in parallel}\ akcja (i).

Wykonanie tej instrukcji polega na wykonaniu dwóch równoległych operacji:

  • przydzielenie procesora do każdego elementu ze zbioru X,
  • jednoczesne wykonanie przez każdy procesor operacji akcja(i).

Przeważnie zapis i \in X jest w rodzaju ``1\leq i\leq n jeśli X jest zbiorem liczb naturalnych.

Podstawowe typy maszyny PRAM

Mamy kilka rodzajów maszyny PRAM w zależności od konfliktów czytania/zapisu we wspólnej pamięci. Litera C (od ang. concurrent) oznacza możliwość jednoczesnego wykonania operacji przez wiele procesorów, E (od ang. exclusive) wyklucza taką możliwość. Operacjami są R (czytanie, od ang. read) oraz W (zapis, od ang. write) w tej samej komórce przez wiele procesorów w tym samym momencie. Mamy zatem EREW PRAM, CREW PRAM, CRCW PRAM (modelu ERCW nie rozważamy jako zupełnie sztuczny).

Podstawowym naszym modelem PRAM-u będzie CREW PRAM: wiele procesrów może jednocześnie czytać z tej samej komórki, ale tylko jeden może zapisywać.


Prostym przykładem obliczenia na CREW jest liczenie kolejnych wierszy trójkąta Pasacala. Początkowo zakładamy, że A = [0,0,0,0,0,1]. Wykonujemy:

   repeat 6 times

   for each 1 \le i \le 5 do in parallel

      A[i]:=A[i]+A[i+1]

Kolejnymi wartościami tablicy A są wektory:

0 0 0 0 0 1
0 0 0 0 0 1
0 0 0 0 1 1
0 0 0 1 2 1
0 0 1 3 3 1
0 1 4 6 4 1
1 5 10 10 5 1

Klasa NC i równoległe algorytmy optymalne

Najważniejszą klasą problemów algorytmicznych stanowią problemy, które można obliczyć w czasie wielomianowo-logarytmicznym, używając wielomianowej liczby procesorów. Klasę tę oznaczamy przez NC, odpowiadające algorytmy nazywamy algorytmami typu NC.

Niech \cal{P} oznacza klasę problemów wykonywanych w deterministycznym czasie wielomianowym na maszynie sekwencyjnej. Podstawowym problemem teoretycznym algorytmiki równoległej jest pytanie:

\cal{P}\ = \ NC \?


Podobnie jak problemy NP-zupełne można zdefiniować problemy P-zupełne. Są to te problemy X \in \cal{P}, takie że dla każdego innego problemu Y \in \cal{P} istnieje NC-redukcja Y do X. Inaczej mówiąc

\cal{P} =  NC \ wtedy i tylko wtedy gdy X \in NC


Przykłady problemów P-zupełnych: programowanie liniowe, maksymalny przepływ w grafie, konstrukcja drzewa DFS, obliczanie wartości układów logicznych, sprawdzanie czy gramatyka bezkontekstowa generuje język pusty.

Przez pracę algorytmu równoległego rozumiemy liczbę procesorów pomnożoną przez czas. Algorytm jest optymalny, gdy jego praca jest tego samego rzędu co czas najlepszego znanego algorytmu sekwencyjnego dla danego problemu.

W szczególności interesują nas algorytmy, które są jednocześnie optymalne i są typu NC. Z praktycznego punku widzenia czynnik log n przy liczbie procesorów i przy pracy algorytmu nie jest zbyt istotny. Natomiast czynnik logarytmiczny jest istotny, jeśli chodzi o równoległy czas. W tym przypadku potęga logarytmu odgrywa podobną rolę co potęga wielomianu opisującego czas obliczenia sekwencyjnego. Różnica między CRCW i CREW PRAMem w aspekcie klasy NC polega przeważnie na dodaniu jednego czynnika logarytmicznego w funkcji równoległego czasu obliczenia.

Przykłady obliczeń na modelu CRCW PRAM

W przypadku CRCW PRAM założymy, że procesory, jeśli wpisują jednocześnie do tej samej komórki pamięci, to wpisują to samo. Na przykład jeśli początkowo output=0, wówczas następujący algorytm obliczy logiczną alternatywę w czasie stałym na CRCW PRAM.


for each 1 \le i \le n do in parallel

    if A[i] =1 then output=1;


Na CREW potrzebujemy logarytmicznego czasu równoległego, aby to zrobić. Pokażemy jeszcze dwa proste problemy, które można na CRCW PRAM wykonać w czasie stałym. Następujący algorytm liczy pierwszą pozycję minimalnego elementu w tablicy C[1 . . n] w czasie O(1).

Algorytm


for each 1 \le i\le n do in parallel
   M[i]:= 0 ;

for each 1 \le i,\ j \le n do in parallel
    if i \ne j and C[i] \le C[j] then M[i]:=1;

for each 1 \le i\le n do in parallel
   if M[i] = 0 then output :=i ;

Oznaczmy ten algorytm przez A_1. Algorytm korzysta z n^2 procesorów. Spróbujemy zmniejszyć tę liczbę do O(n^{1+\epsilon}), zachowując czas O(1) dla dowolnie małego \epsilon >0.

Niech

P_k(n)= n^{1+\epsilon_k}, gdzie \epsilon_k\ =\ \frac{1}{2^{k}-1}.


Przypuśćmy, że mamy algorytm A_k liczenia minimum w czasie O(1) z O(P_k(n)) procesorami. Skonstruujemy algorytm A_{k+1}, który działa w czasie stałym i używa tylko O(P_{k+1}(n)) procesorów.


Algorytm A_{k+1}



niech \alpha=\frac{1}{2^{k}+1};

podziel tablicę C na rozłączne bloki rozmiaru n^{\alpha} każdy;

równolegle zastosuj algorytm A_k do każdego z tych bloków;

zastosuj algorytm A_k do tablicy C' składającej się z \frac{n}{n^{\alpha}} minimów w blokach.


Algorytm A_{k+1} działa w czasie O(1), korzystając z P_{k+1}(n)) procesorów. Uzasadnienie pozostawiamy jako ćwiczenie.

Algorytmy A_2,\ A_3,\  A_4, \ldots używają odpowiednio następującą (asymptotycznie) liczbę procesorów

n^{1+\frac{1}{3}},\ \ n^{1+\frac{1}{15}},\  \ n^{1+\frac{1}{255}},  \ \ n^{1+\frac{1}{65535}} ...

Rozważmy jeszcze na CRCW PRAM następujący problem pierwszej jedynki: dana jest tablica zerojedynkowa, znaleźć pozycję pierwszej jedynki (od lewej).

Następujący algorytm rozwiązuje problem w czasie stałym z kwadratową liczbą procesorów. Zakładamy na razie, że w ciągu jest jakaś jedynka.

Algorytm Pierwsza-Jedynka-1



for each 1 \leq i < j \le n do in parallel

   if A[i]=1 and A[j]=1 then A[j]:=0;

for each 1 \leq i \le n do in parallel

   if A[i]=1 then FirstOne :=i.


Możemy podobnie łatwo sprawdzić, czy w ogóle jest jedynka.

Algorytm CzyJestJedynka



jest-jedynka :=0;

for each 1 \leq i \le n do in parallel

   if A[i]=1 then jest-jedynka :=1;


Oba powyższe algorytmy korzystają z O(n^2) procesorów. Możemy łatwo tę liczbę zmniejszyć do liniowej.


Algorytm Pierwsza-Jedynka


(1) Podziel tablicę A na segmenty długośći \sqrt{n};

(2) W każdym segmencie zastosuj algorytm CzyJestJedynka;

(3) Otrzymujemy ciąg zerojedynkowy C długości \sqrt{n} jako wynik kroku (2);

(4) znajdź pierwszą jedynkę w ciągu C za pomocą algorytmu Pierwsza-Jedynka-1;

(5) Zastosuj algorytm Pierwsza-Jedynka-1 do segmentu odpowiadającego pierwszej jedynce w C;

W ten sposób stosujemy trzy razy algorytm o pracy kwadratowej do segmentu długości \sqrt{n}. Otrzymujemy złożoność O(\sqrt{n}^2) = O(n). Czas jest O(1).

Drzewa i równoległa wersja metody dziel i zwyciężaj

Do szybkich obliczeń równoległych najbardziej nadają się problemy związane z drzewami, chociaż czasami w tych problemach nie widać od razu struktury drzewiastej. Struktura taka odpowiada również drzewu rekursji. Jako przykład rozważmy problem obliczenia sumy A[1]+A[2]+\ldots A[n]. Dla uproszczenia załóżmy, że n jest potęgą dwójki.

Grafika:Parallel1-2.png
Rysunek 2: Metoda pełnego zrównoważonego drzewa binarnego: układ arytmetyczny obliczania sumy. Maksymalny poziom m=\log n.

Wysokością węzła jest jego maksymalna odległość od liścia, wysokość liścia wynosi 0. Przez p-ty poziom rozumiemy zbiór węzłów o wysokości p. Załóżmy, że elementy A[1], A[2], ..A[n] są umieszczone w liściach pełnego zrównoważonego drzewa binarnego, następnie wykonujemy (patrz rysunek):


for p:=1 to \log n do

    oblicz jednocześnie wartości węzłów na poziomie p-tym;


Drzewo jest strukturą koncepcyjną, każdemu węzłowi możemy przypisać miejsce w pamięci. W naszym przypadku węzły na poziomie p-tym mogą odpowiadać elementom

A[2^p],\ A[2*2^p], .. A[3*2^p]\ldots.

Poprzedni algorytm można zapisać w formie:

for p:=1 to \log n do

   \Delta=2^p;

   for each 1 \le i \le n/\Delta do in parallel

       A[i*\Delta]\ :=\  A[i*\Delta]-\Delta/2]+A[i*\Delta];

wynik := A[n];


Grafika:Parallel1-3.png
Rysunek 3. Koncepcyjna struktura równoległej wersji metody dziel i zwyciężaj.


Drzewa odpowiadają w pewnym sensie rekursji. Wysokość drzew odpowiada czasowi równoległemu. Podobnie głębokość rekursji odpowiada czasowi równoległemu. Równoległa wersja tzw. metody dziel i zwyciężaj polega na tym, że wywołania rekurencyjne wykonujemy jednocześnie (patrz rysunek). Oczywiście może się zdarzyć, że jedno z nich zakończy się wcześniej. Niemniej algorytm czeka na zakończenie obu wywołań.


Pokażemy dwa przykłady zastosowania metody dziel i zwyciężaj w wersji równoległej. Zaczniemy od sumy elementów tablicy A[1..n]. Wynik obliczamy jako SUMA(1,n). Zakładamy znowu, że n jest potęgą dwójki.

funkcja SUMA(1,n)

if j=i then return A[i] else

   do in parallel

      wynik1 := SUMA(i,\ \lfloor (i+j)/2 \rfloor;

      wynik2 := SUMA(\lceil (i+j)/2 \rceil,\ j);

   return wynik1 + wynik2;


Podobny jest schemat sortowania na PRAMie. Niech ParallelMerge(x) będzie algorytmem, który, otrzymawszy tablicę x z posortowanymi lewą i prawą połową, da w wyniku posortowaną tablicę x. Łatwo to zrobić w czasie O(\log n) z n procesorami. Dla i>n/2-ty procesor znajduje w pierwszej połówce sekwencyjnie, metodą binary search, najmniejszy element większy od x[i]. Wymaga to czasu O(\log n). Potem każdy procesor{\em wie} gdzie wstawić {\em swój} element. W sumie otrzymujemy algorytm sortowania w czasie O(\log^2) z n procesorami.


funkcja ParallelSort(x);

n:=size(x);

if n>1 then

    do in parallel

       ParallelSort(FirstHalf(x));

       ParallelSort(SecondHalf(x));

   ParallelMerge(x)


Liczbę procesorów można zmniejszyć do n/(\log n)). Natomiast nietrywialnym jest zmniejszenie czasu na CREW PRAM. Zostało to zrobione przez Richarda Cole'a, który skonstruował algorytm działający w czasie O(\log n) z O(n) procesorami maszyny EREW PRAM. Algorytm ten jest bardzo interesujący, ale niestety również bardzo skomplikowany.

Układy arytmetyczne: sumy prefiksowe

Być może najbardziej podstawowym modelem obliczeń równoległych są układy arytmetyczne (lub logiczne): acykliczne grafy z przypisaniem pewnych operacji węzłom wewnętrznym. Każdy węzeł liczy pewną wartość w momencie, gdy wartości jego poprzedników są policzone. Podobnie jak w drzewie, możemy zdefiniować pojęcie liścia: węzeł bez poprzedników. Natomiast graf nie musi mieć jednego korzenia. Zamiast korzenia w grafie wyróżniamy węzły wynikowe (na rysunku te, z których wychodzi strzałka do ‘’nikąd’’).

Równoległy czas obliczenia odpowiada maksymalnej wysokości węzła. Poziomy definiujemy podobnie jak dla drzewa. Algorytm równoległy w jednym równoległym kroku oblicza kolejny poziom. Liczba procesorów odpowiada z reguły maksymalnemu rozmiarowi poziomu, chociaż możemy inaczej rozplanować obliczenie, gdy jedne poziomy są duże, a drugie małe (ale trzeba wtedy zmienić strukturę grafu tak, aby temu odpowiadała).

Przykładem układu arytmetycznego jest drzewo z rysunku powyżej, które opisuje sumowanie n elementów.

Zajmiemy się teraz pewnym rozszerzeniem problemu sumowania. Niech \oplus będzie pewną łączną operacją arytmetyczną (np. suma, mnożenie, maksimum, minimum, pozycja pierwszej jedynki z lewej strony, podobnie z prawej strony).

Problem sum prefiksowych. dany wektor x rozmiaru n, obliczyć ciąg y taki, gdzie

y[1]=x[1], y[2]=x[1]\oplus x[2], y[3]=x[1]\oplus x[2]\oplus x[3], ...

gdzie


Opiszemy dwa rekurencyjne algorytmy dla tego problemu. Niech FirstHalf, SecondHalf oznaczają lewą i prawą (odpowiednio) połówkę ciągu. Zakładamy, że n jest potęgą dwójki.

Algorytm PrefSums1(x);



n:=size(x);

ifn>1 then

   do in parallel

      PrefSums1(FirstHalf(x));

      PrefSums1(SecondHalf(x));

   for each n/2 < j\leq n, do in parallel

      x[j]:= x[n/2] \oplus x[j];

Układ arytmetyczny odpowiadający algorytmowi PrefSums1 jest przedstawiony na Rysunku 4 dla n=4 i n=8. Zauważmy, że zasadniczą częścią układu dla n=8 są dwie kopie układu dla n=4. Dodatkowo dodajemy węzły odpowiadającej ostatniej instrukcji w algorytmie PrefSums1.


Grafika:Parallel1-4.png
Rysunek 4. Układ arytmetyczny odpowiadający algorytmowi PrefSums1. Kolejne grafy powstają jako podwójne kopie porzednich grafów (dla n/2 elementów) z dodanymi elementami odpowiadającymi operacji x[j]:= x[n/2] \oplus x[j].


Opiszemy teraz inny algorytm rekurencyjny, w którym mamy tylko jedno wywołanie rekurecyjne (w jednej instancji rekursji).

Algorytm PrefSums2(x)



n:=size(x);

if n>1 then

   utwórz nową tablicę y;

   for each 1 \le i \le n/2 do in parallel

       y[i]\ :=\ x[2i-1]\oplus x[2i];

   PrefSums2(y);

   for each 1 \le i \le n/2 do in parallel

      x[2i]:= y[i];       if i>1 then x[2i-1]:= y[i-1] \oplus x[2i-1];


Układ arytmetyczny odpowiadający algorytmowi jest pokazany na rysunku 5.


Grafika:Parallel1-5.png
Rysunek 5. Układ arytmetyczny odpowiadający PrefSums2. Kolejny graf składa się z pojedynczej kopii poprzedniego grafu (dla n/2) oraz n/2-1 dodatkowych operacji.