Zaawansowane algorytmy i struktury danych/Wykład 13

Algorytmy równoległe I

W module tym zajmiemy się przyspieszaniem obliczeń za pomocą korzystania z wielu procesró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 niskopoziomowym, ale niewątpliwie bardzo istotnym.

Model równoległej abstrakcyjnej Maszyny PRAM

Na początku rozważymy wyidealizowany model obleczeń 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).

Rysunek 1: Struktura koncepycyjna 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 komputerm 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ówoległość jest wyrażona poprzez następującą instrukcję:

$\text{<mrow data-mjx-texclass="ORD">forall</mrow>}<math>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\le i\le n$ jeśli $X$ jest zbiorem liczb naturalnych.

Podstawowe typy maszyny PRAM

Mamy kilka rodziajów maszyny PRAM w zależności od konfliktów czytania/zapisu we wspólnej pamięi. 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 mdoelem PRAMu będzie CREW PRAM: wiele procesrów może jednocześnie czytać z tej samej komórki, ale tylko jeden może zapisywać. \myskip 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: \vskip 0.2cm \begin{center} \begin{minipage}{10cm} repeat 6 times

for each $1\le i\le 5$ do in parallel \\ \hspace*{1cm} $A[i]:=A[i]+A[i+1]$ \end{minipage} \end{center} \vskip 0.2cm \noindent Kolejnymi wartościami tablicy A są wektory: \vskip 0.2cm \begin{center} \noindent 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 1

\noindent 0\ \ 0\ \ 0\ \ 0\ \ 0\ \ 1

\noindent 0\ \ 0\ \ 0\ \ 0\ \ 1\ \ 1

\noindent 0\ \ 0\ \ 0\ \ 1\ \ 2\ \ 1

\noindent 0\ \ 0\ \ 1\ \ 3\ \ 3\ \ 1

\noindent 0\ \ 1\ \ 4\ \ 6 \ \ 4\ \ 1

\noindent 1\ \ 5\ 10 \ 10\ \ 5\ \ 1 \end{center}

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 ${𝒫}$ oznacza klasę problemów wykonywanych w deterministycznym czasie wielomianowym na maszynie sekwencyjnej. Podstawowym problemem teoretycznym algorytmiki równoległej jest pytanie:

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

\centerline{${𝒫}{=}{𝒩}{𝒞}$ wtedy i tylko wtedy gdy \ $X\in NC$} \myskip 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. \myskip Przez pracę algorytmu równoległego rozumiemy liczbę procesorów pomnożoną orzez 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 $logn$ 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. 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. \myskip for each $1\le i\le n$ do in parallel \\ \hspace*{1cm} if $A[i]=1$ then output=1; \vskip 0.6cm 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 czasei stałym. Następujący algorytm liczy pierwszą pozycję minimalnego elementu w tablicy C[1 . . n] w czasie O(1). \vskip 0.4cm for each $1\le i\le n$ do in parallel $M[i]$ := 0;

for each $1\le i,j\le n$ do in parallel \\ \hspace*{1cm} if $i\ne j$ and $C[i]\le C[j]$ then $M[i]$:=1;

for each $1\le i\le n$ do in parallel \\ \hspace*{1cm} if $M[i]$ = 0 then output :=i ; \myskip Oznaczmy ten algorytm przez $A_1$. Algorytm korzysta z $n^2$ procesorów. Spróbujemy zmniejszyć tę liczbę do $O(n^{1+ϵ})$ zachowując czas $O(1)$, dla dowolnie małego $ϵ>0$.

\noindent Niech

\centerline{$P_k(n)=n^{1+{ϵ}_k}$, gdzie ${ϵ}_k=\frac{1}{2^k-1}$.} \myskip \noindent 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. \vskip 0.5cm \noindent Algorytm}\ Parser nie mógł rozpoznać (błąd składni): {\displaystyle 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. \vskip 0.5cm \noindent Algorithm $A_{k+1}$ działa w czasie $O(1)$ korzystając z $P_{k+1}(n))$ procesorów. Uzasadnienie pozostawiamy jako ćwiczenie. \vskip 0.3cm \noindent Algorytmy $A_2,A_3,A_4,\dots$ 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}}..$. \myskip Rozważmy jeszcze na CRCW PRAM następujący problem pierwszej jedynki: \ dana tablica zerojedynkowa, znaleźć pozycję pierwszej jedynki (od lewej). \myskip 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ąs jedynka. \vskip 0.3cm Algorytm Pierwsza-Jedynka-1;

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

\hspace*{0.6cm} if A[i]=1] and A[j]=1 then A[j]:= 0;

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

\hspace*{0.6cm} if A[i]=1 then FirstOne :=i. \vskip 0.3cm \noindent Możemy podobnie łatwo sprawdzić czy w ogóle jest jedynka. \vskip 0.3cm Algorytm CzyJestJedynka;

jest-jedynka := 0;

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

\hspace*{0.6cm} if A[i]=1 then jest-jedynka := 1;

\vskip 0.2cm \noindent Oba powyższe algorytmy korzystają z $O(n^2)$ procesorów. Możemy łatwo tę liczbę zmniejszyć do liniowej. \myskip Algorytm Pierwsza-Jedynka; \vskip 0.1cm (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 Parser nie mógł rozpoznać (błąd składni): {\displaystyle \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 \\ \hspace*{1.4cm} pierwszej jedynce w C; \myskip 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)$. 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]+\dots A[n]$. Dla uproszczenia załóżmy, że n jest potęgą dwójki.

\begin{figure}[htbp] \begin{center} \mbox{\ } \includegraphics[width=10.cm]{parallel_fig2.eps} \caption{ Metoda pełnego zrównoważonego drzewa binarnego:\ układ arytmetyczny obliczania sumy. Maksymalny poziom $m=\log n$.}

\end{center} \end{figure}

\noindent 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 drrzewa binarnego, następnie wykonujemy (patrz rysunek): \vskip 0.4cm for $p:=1$ to $\log n$ do

\hspace*{0.5cm} oblicz jedmocześnie wartości węzłów na poziomie $p$-tym; \myskip 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 \\ \centerline{$A[2^p],A[2*2^p],..A[3*2^p]\dots$. } Poprzedni algorytm można zapisać w formie: \vskip 0.4cm \noindent for $p:=1$ to $\log n$ do\\ \hspace*{0.5cm} $\mathrm{\Delta }=2^p$; \\ \hspace*{0.5cm} for each $1\le i\le n/\mathrm{\Delta }$ do in parallel \\ \hspace*{0.9cm} $A[i*\mathrm{\Delta }]:=A[i*\mathrm{\Delta }]-\mathrm{\Delta }/2]+A[i*\mathrm{\Delta }]$; \\ wynik := $A[n]$; \myskip \begin{figure}[hbtp] \begin{center} \mbox{\ } \includegraphics[width=9.cm]{parallel_fig3.eps} \caption{Koncepcyjna struktura równoległej wersji metody {\em dziel i zwycieżaj}. }

\end{center} \end{figure}

Drzewa odpowiadają w pewnym sensie rekursji. Wysokość drzew odpowiada czasowi równoległemu. Podobnie głębokość rekursji odpowiada też 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. Tym 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. \myskip \begin{center} \begin{minipage}{12cm} \vskip0.3cm \hspace*{0.6cm}\textbf{funkcja} $SUMA(i,j)$; \\ \hspace*{1.2cm}\textbf{if} $j=i$ \textbf{then return} $A[i]$ \textbf{else }\\ \hspace*{1.8cm}{do in parallel} \\ \hspace*{2.5cm} wynik1 := $SUMA(i,\lfloor (i+j)/2\rfloor$;\\ \hspace*{2.5cm} wynik2 := $SUMA(\lceil (i+j)/2\rceil ,j)$;\\ \hspace*{1.8cm}\textbf{return} wynik1 + wynik2; \vskip0.4cm \end{minipage} \end{center} \myskip \noindent 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 wyniku tablicę $x$ posortowaną. łatwo to zrobić w czasie $O(\log n)$ z $n$ procesorami. Dla $i>n/2$-ty procesor znajduje w pierwszej połówce sekwencyjnie metodą {\em 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. \begin{center} \begin{minipage}{12cm} \vskip0.3cm \hspace*{0.6cm}\textbf{funkcja} $ParallelSort(x)$; \\ \hspace*{1.2cm}$n:=size(x)$;\\ \hspace*{1.2cm}\textbf{if} $n>1$ \textbf{then } \\ \hspace*{1.8cm}\textbf{do in parallel}\\ \hspace*{2.5cm} $ParallelSort(FirstHalf(x))$;\\ \hspace*{2.5cm} $ParallelSort(SecondHalf(x))$;\\ \hspace*{1.8cm} $ParallelMerge(x)$ \vskip0.4cm \end{minipage} \end{center} \myskip Liczbę procesorów można zmniejszyc 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 czaie $O(\log n)$ z $O(n)$ procesorami maszyny EREW PRAM. Algorytm ten jest bardzo interesujący ale skomplikowany. \newpage Być może najbardziej podstawowym modelem obliczeń równoległych są układy arytmetyczne (lub logiczne): acyckliczne 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 zdefinować 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 {\em do nikąd}).

Równoległy czas obliczenia odpowiada maksymalej 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 roziarowi poziomu, chociaż możemy inaczej rozplanować obliczenie gdy jedne poziomy są duże, a drugie małe (ale możemy 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).

\noindentProblem sum p[refiksowych.\ dany wektor $x$ rozmiaru $n$, obliczyć ciąg $y$ taki, gdzie \myskip \begin{center} $y[1]=x[1]$, $y[2]=x[1]\oplus x[2]$, $y[3]=x[1]\oplus x[2]\oplus x[3]$, \ldots \end{center} gdzie \myskip Opiszemy dwa rekurencyjne algorytmy dla tego problemu. Niech $FirstHalf$, $SecondHalf$ oznaczją lewą i prawą (odpowiednio) połówkę ciągu. Zakładamy, że n jest potęgą dwójki. \myskip \begin{center} \begin{minipage}{12cm} \vskip0.3cm \hspace*{0.6cm}\textbf{Algorytm} $PrefSums1(x)$; \\ \hspace*{1.2cm}$n:=size(x)$;\\ \hspace*{1.2cm}\textbf{if} $n>1$ \textbf{then } \\ \hspace*{1.8cm}\textbf{do in parallel}\\ \hspace*{2.5cm} $PrefSums1(FirstHalf(x))$;\\ \hspace*{2.5cm} $PrefSums1(SecondHalf(x))$;\\ \hspace*{1.8cm}\textbf{for each } $n/2<j\le n$, \textbf{do in parallel} \\ \hspace*{2.4cm} $x[j]:=x[n/2]\oplus x[j]$;\\ \vskip0.4cm \end{minipage} \end{center} \myskip

Układ arytmetyczny odpowiadający algorytmowi PrefSums1 jest przedstawiony na Rysunku #parallel_fig5 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.

\begin{figure}[bhtp] \begin{center} \mbox{\ } \includegraphics[width=11.5cm]{parallel_fig5.eps} \caption{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]$. }

\end{center} \end{figure}

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

\myskip \begin{center} \begin{minipage}{12cm} \vskip0.3cm \hspace*{0.6cm}\textbf{Algorytm} $PrefSums2(x)$; \\ \hspace*{1.2cm}$n:=size(x)$;\\ \hspace*{1.2cm}\textbf{if} $n>1$ \textbf{then } \\ \hspace*{1.8cm} utwórz nową tablicę y;\\ \hspace*{1.8cm} for each } $1\le i\le n/2$ \textbf{do in parallel\\ \hspace*{2.5cm} $y[i]:=x[2i-1]\oplus x[2i]$;\\ \hspace*{1.8cm} $PrefSums2(y)$;\\ \hspace*{1.8cm} for each } $1\le i\le n/2$ \textbf{do in parallel\\ \hspace*{2.4cm} $x[2i]:=y[i]$;\\ \hspace*{2.4cm} if $i>1$ then \ $x[2i-1]:=y[i-1]\oplus x[2i-1]$;\\ \vskip0.4cm \end{minipage} \end{center} \myskip Układ arytmetyczny odpowiadający algorytmowi jest pokazany na rysunku  #parallel_fig4. \begin{figure}[bhtp] \begin{center} \mbox{\ } \includegraphics[width=7.2cm]{parallel_fig4.eps} \caption{Układ arytmetyczny odpowiadający PrefSums2. Kolejny graf składa się z pojedyńczej kopii poprzedniego grafu (dla $n/2$), oraz $n/2-1$ dodatkowych operacji. }

\end{center} \end{figure}