|
|
| Linia 1: |
Linia 1: |
|
| |
|
| ''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki''
| |
|
| |
| ==Rozwiązywanie układów równań liniowych==
| |
|
| |
| Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń
| |
| zmiennoprzecinkowych, przy czym wiele znanych, ''matematycznie równoważnych'' metod
| |
| rozwiązywania takich zadań, ma ''diametralnie różne'' własności numeryczne.
| |
| Bardzo ważną klasą takich zadań jest rozwiązywanie układów równań liniowych
| |
|
| |
| <center><math>\displaystyle
| |
| Ax = b,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle A</math> jest nieosobliwą macierzą <math>\displaystyle N\times N</math>, a dany wektor prawej strony <math>\displaystyle b\in R^N</math>.
| |
|
| |
| W
| |
| praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
| |
| specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
| |
|
| |
| Rozwiązywanie układów równań liniowych jest sercem wielu innych algorytmów
| |
| numerycznych, dlatego nie dziwi, że szacuje się, że około 75 procent czasu
| |
| obliczeniowego superkomputerów jest wykorzystywanych właśnie na rozwiązywanie
| |
| takich zadań.
| |
|
| |
| Okazuje się, że kilka znanych w matematyce sposobów rozwiązywania układów równań
| |
| liniowych, takie jak:
| |
| * metoda wyznacznikowa (wzory Cramera)
| |
| * obliczenie macierzy <math>\displaystyle A^{-1}</math> i następnie <math>\displaystyle x = A^{-1}b</math>
| |
|
| |
| '''nie nadają się''' do numerycznego rozwiązywania takich zadań.
| |
|
| |
| O tym, jak ''skutecznie'' rozwiązywać takie zadania, jakie jest ich
| |
| uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z
| |
| dobrych metod jest metoda eliminacji Gaussa.
| |
|
| |
| ===Proste układy równań===
| |
|
| |
| Niektóre układy równań można bardzo łatwo rozwiązać. Zgodnie z zasadą
| |
|
| |
| <blockquote>
| |
| Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
| |
| </blockquote>
| |
|
| |
| w dalszej kolejności pokażemy, jak dowolny układ równań sprowadzić do sekwencji
| |
| dwóch (czasem trzech) łatwych do rozwiązania układów równań. Ale... jakie układy
| |
| równań są "łatwe"?
| |
|
| |
| ====Układy z macierzą trójkątną====
| |
|
| |
| Rozważmy układ z macierzą
| |
| trójkątną <math>\displaystyle A</math>. Będą nas szczególnie interesować macierze
| |
| ''trójkątne górne'', dla których <math>\displaystyle a_{i,j}=0</math> gdy <math>\displaystyle i>j</math>, oraz
| |
| macierze ''trójkątne dolne'' z jedynkami na przekątnej, tzn.
| |
| <math>\displaystyle a_{i,j}=0</math>, <math>\displaystyle i<j</math>, oraz <math>\displaystyle a_{i,i}=1</math>. Macierze pierwszego rodzaju
| |
| będziemy oznaczać przez <math>\displaystyle U</math>, a drugiego rodzaju przez <math>\displaystyle L</math>.
| |
|
| |
| <center><math>\displaystyle L = \begin{pmatrix}
| |
| 1 & & & & & \\
| |
| * & 1 & & & & \\
| |
| * & * & 1 & & & \\
| |
| * & * & * & 1 & & \\
| |
| \vdots & \vdots & \vdots & \ddots & \ddots & \\
| |
| * & * & * & \cdots & * & 1
| |
| \end{pmatrix} ,
| |
| \qquad
| |
| U = \begin{pmatrix}
| |
| * & * & * & * & \cdots & * \\
| |
| & * & * & * & \cdots & * \\
| |
| & & * & * & \cdots & * \\
| |
| & & & * & \ddots & \vdots \\
| |
| & & & & \ddots & * \\
| |
| & & & & & * \end{pmatrix}
| |
| </math></center>
| |
|
| |
| Układ z nieosobliwą macierzą trójkątną górną
| |
|
| |
| <center><math>\displaystyle
| |
| U\, x\;=\; c,
| |
| </math></center>
| |
|
| |
| <math>\displaystyle U=(u_{i,j})</math>, <math>\displaystyle c=(c_j)</math>, można rozwiązać stosując algorytm:
| |
|
| |
| {{algorytm|Podstawienie w tył||
| |
| <pre>
| |
|
| |
| <math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
| |
| for (i <nowiki>=</nowiki> N-1; i ><nowiki>=</nowiki> 1; i--)
| |
| <math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
| |
| u_{i,j}x_j^*\right)/u_{i,i}</math>;
| |
| </pre>}}
| |
|
| |
| (Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy
| |
| implikuje, że <math>\displaystyle u_{i,i}\ne 0</math>, <math>\displaystyle \forall i</math>.) Podobnie, układ
| |
| <math>\displaystyle L x= c</math> rozwiązujemy algorytmem:
| |
|
| |
| {{algorytm|Podstawienie w przód||
| |
| <pre>
| |
|
| |
| <math>\displaystyle x_1 = c_1</math>;
| |
| for (i<nowiki>=</nowiki>2; i <<nowiki>=</nowiki> N; i++)
| |
| <math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;
| |
| </pre>}}
| |
|
| |
| Oba algorytmy wymagają rzędu <math>\displaystyle N^2/2</math> mnożeń lub dzieleń i
| |
| <math>\displaystyle N^2/2</math> dodawań lub odejmowań, a więc łącznie <math>\displaystyle O(N^2)</math>
| |
| działań arytmetycznych.
| |
|
| |
| ====Układy z macierzą ortogonalną====
| |
|
| |
| Równie tanio można rozwiązać układ równań
| |
|
| |
| <center><math>\displaystyle
| |
| Q x = b,
| |
| </math></center>
| |
|
| |
| gdy <math>\displaystyle Q</math> jest macierzą ortogonalną, to znaczy <math>\displaystyle Q^TQ = I</math>. Rzeczywiście, z
| |
| ortogonalności mamy natychmiast, że
| |
|
| |
| <center><math>\displaystyle
| |
| x = Q^T b
| |
| </math></center>
| |
|
| |
| i w konsekwencji <math>\displaystyle x</math> można wyznaczyć kosztem takim jak koszt mnożenia macierzy
| |
| przez wektor, czyli <math>\displaystyle O(N^2)</math> operacji.
| |
|
| |
| Podobnie, gdy <math>\displaystyle Q\in C^{N\times N}</math> jest unitarna, to znaczy <math>\displaystyle Q^HQ = I</math>,
| |
| rozwiązaniem układu równań jest
| |
|
| |
| <center><math>\displaystyle
| |
| x = Q^H b.
| |
| </math></center>
| |
|
| |
| ===Metoda eliminacji Gaussa===
| |
|
| |
| W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym
| |
| rozwiązywania układu równań
| |
|
| |
| <center><math>\displaystyle Ax=b</math></center>
| |
|
| |
| okazuje się popularna
| |
| eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm
| |
| ten wyrazimy w języku tzw. ''rozkładu LU'' macierzy, to znaczy,
| |
| sprowadzającego zadanie do znalezienia
| |
| macierzy trójkątnej dolnej <math>\displaystyle L</math> (z jedynkami na diagonali) oraz trójkątnej górnej
| |
| <math>\displaystyle U</math> takich, że
| |
| <center><math>\displaystyle
| |
| A = LU,
| |
| </math></center>
| |
|
| |
| a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:
| |
|
| |
| {{algorytm|Rozwiązanie układu z wykorzystaniem rozkładu LU||
| |
| <pre>
| |
|
| |
| Znajdź rozkład <math>\displaystyle A=LU</math>;
| |
| Rozwiąż <math>\displaystyle Ly = b</math> przez podstawienie w przód;
| |
| Rozwiąż <math>\displaystyle Ux = y</math> przez podstawienie w tył;
| |
| </pre>}}
| |
|
| |
| Przypuśćmy, że taki rozkład <math>\displaystyle A=LU</math> istnieje. Wówczas, zapisując macierze w postaci
| |
| blokowej, eskponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy,
| |
| mamy
| |
|
| |
| <center><math>\displaystyle
| |
| \begin{pmatrix}
| |
| a_{11} & a_{12}^T\\
| |
| a_{21} & A_{22}
| |
| \end{pmatrix}
| |
| =
| |
| \begin{pmatrix}
| |
| 1 & 0^T\\
| |
| l_{21} & L_{22}
| |
| \end{pmatrix}
| |
| \begin{pmatrix}
| |
| u_{11} & u_{12}^T\\
| |
| 0 & U_{22},
| |
| \end{pmatrix}
| |
| </math></center>
| |
|
| |
| skąd (mnożąc blokowo macierz <math>\displaystyle L</math> przez <math>\displaystyle U</math>) wynika, że
| |
| * <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest
| |
| kopią pierwszego wiersza <math>\displaystyle A</math>,
| |
| * <math>\displaystyle l_{21} = a_{21}/u_{11}</math>, więc pierwsza kolumna <math>\displaystyle L</math> powstaje przez
| |
| podzielenie wszytkich elementów wektora <math>\displaystyle a_{21}</math> przez element na diagonali
| |
| <math>\displaystyle a_{11}</math>,
| |
| * <math>\displaystyle A_{22} - l_{21}u_{12}^T = L_{22}U_{22}</math>, a więc znalezienie podmacierzy
| |
| <math>\displaystyle L_{22}</math> oraz <math>\displaystyle U_{22}</math> sprowadza się do znalezienia rozkładu LU zmodyfikowanego
| |
| bloku <math>\displaystyle A_{22}</math> macierzy <math>\displaystyle A</math>,
| |
| wymiaru <math>\displaystyle (N-1)\times (N-1)</math>.
| |
|
| |
| Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie
| |
| rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
| |
| klasycznych pętli. Jest to ważne w praktyce numerycznej, gdyż rekurencja
| |
| kosztuje: zarówno pamięć, jak i czas.
| |
|
| |
| Ponadto zauważmy, że opisany algorytm możemy wykonać w miejscu, nadpisując
| |
| elementy <math>\displaystyle A</math> elementami macierzy <math>\displaystyle U</math> i <math>\displaystyle L</math> (jedynek z diagonali <math>\displaystyle L</math> nie musimy
| |
| pamiętać, bo wiemy a priori, że tam są).
| |
|
| |
| {{algorytm|Rozkład LU metodą eliminacji Gaussa||
| |
| <pre>
| |
|
| |
| for k<nowiki>=</nowiki>1:N-1
| |
| if <math>\displaystyle a_{kk}</math> <nowiki>=</nowiki><nowiki>=</nowiki> 0
| |
| STOP;
| |
| end
| |
| for i<nowiki>=</nowiki>k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| |
| <math>\displaystyle a_{ik}</math> <nowiki>=</nowiki> <math>\displaystyle a_{ik}/a_{ii}</math>;
| |
| end
| |
| for j<nowiki>=</nowiki>k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
| |
| for i<nowiki>=</nowiki>k+1:N
| |
| <math>\displaystyle a_{ij}</math> -<nowiki>=</nowiki> <math>\displaystyle a_{ik}a_{kj}</math>;
| |
| end
| |
| end
| |
| end
| |
| </pre>}}
| |
|
| |
| Łatwo przekonać się, że <math>\displaystyle k</math>-ty obrót zewnętrznej pętli (tzn. <math>\displaystyle k</math>-ty krok
| |
| algorytmu rozkładu LU) kosztuje rzędu <math>\displaystyle 2(N-k)^2</math> operacji arytmetycznych, skąd łączny koszt tego
| |
| algorytmu rozkładu LU wynosi około <math>\displaystyle \frac{4}{3}N^3</math>.
| |
|
| |
| Jeśli więc wykorzystać rozkład LU do rozwiązywania układu równań <math>\displaystyle Ax=b</math>, to mamy
| |
| następujące zestawienie kosztów:
| |
| * Koszt znalezienia rozkładu <math>\displaystyle A=LU</math>: <math>\displaystyle O(N^3)</math>;
| |
| * Koszt rozwiązania układu <math>\displaystyle Ly=b</math>: <math>\displaystyle O(N^2)</math>;
| |
| * Koszt rozwiązania układu <math>\displaystyle Ux=y</math>: <math>\displaystyle O(N^2)</math>.
| |
|
| |
| Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania
| |
| wynosi już tylko <math>\displaystyle O(N^2)</math>.
| |
|
| |
| {{uwaga|Złożoność obliczeniowa zadania rozwiązania układu równań liniowych||
| |
|
| |
| Z powyższego wynika, że łączny koszt rozwiązania równania liniowego wynosi
| |
| <math>\displaystyle O(N^3)</math>. Można zastanawiać się, jaka jest ''najmniejsza możliwa'' liczba
| |
| operacji zmiennoprzecinkowych potrzebnych do rozwiązania układu równań
| |
| liniowych.
| |
|
| |
| Można pokazać \cite{Cormen}, że minimalny koszt rozwiązania układu <math>\displaystyle N</math> równań
| |
| liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch
| |
| macierzy <math>\displaystyle N\times N</math>. Znany jest całkiem prosty algorytm rekurencyjny,
| |
| wyznaczający iloczyn dwóch macierzy kosztem <math>\displaystyle 4.7\cdot N^{log_27} \approx 4.7
| |
| \cdot N^{2.807}</math> (algorytm Strassena). Bardziej skomplikowany (i praktycznie
| |
| nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt
| |
| <math>\displaystyle O(N^{2.376})</math>. Tak więc równania liniowe daje się (teoretycznie) rozwiązać
| |
| kosztem <math>\displaystyle O(N^{2.376})</math>.
| |
|
| |
| Jednak w praktyce nawet prosty algorytm Strassena zazwyczaj nie jest stosowany.
| |
| Wynika to stąd, że ma trochę gorsze własności numeryczne oraz wymaga dużo
| |
| dodatkowej pamięci na przechowywanie pośrednich wyników.
| |
| }}
| |
|
| |
| ===Wybór elementu głównego===
| |
|
| |
| Opisany powyżej algorytm rozkładu LU niestety czasem może się załamać: gdy
| |
| napotka w czasie działania zerowy element na diagonali zmodyfikowanej
| |
| podmacierzy, np. chociaż macierz
| |
|
| |
| <center><math>\displaystyle
| |
| A = \begin{pmatrix} 0 & 1\\ 1 & 0
| |
| \end{pmatrix}
| |
| </math></center>
| |
|
| |
| jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od
| |
| razu zetknie się z dzieleniem przez <math>\displaystyle a_{11}=0</math>. Ale wystarczy zamienić ze sobą
| |
| kolejnością wiersze macierzy <math>\displaystyle A</math> (to znaczy, w układzie równań, zamienić ze sobą
| |
| miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez
| |
| problemu.
| |
|
| |
| W praktyce obliczeniowej, aby uzyskać algorytm o [[|Uzupełnij: możliwie dobrych własnościach
| |
| numerycznych]], wykorzystujemy tzw. strategię ''wyboru elementu głównego w
| |
| kolumnie''.
| |
| Polega to na tym, że zanim wykonamy <math>\displaystyle k</math>-ty krok algorytmu rozkładu LU,
| |
| * szukamy w pierwszej kolumnie podmacierzy <math>\displaystyle A(k:N,k:N)</math> szukamy elementu o
| |
| największym module (taki element, na mocy założenia nieosobliwości macierzy,
| |
| jest niezerowy) --- to jest właśnie element główny
| |
| * zamieniamy ze sobą wiersz pierwszy <math>\displaystyle A(k:N,k:N)</math> z wierszem, w którym
| |
| znajduje się element główny
| |
| * zapamiętujemy dokonaną permutację, bo potem --- gdy już przyjdzie do
| |
| rozwiązywania układu równań --- będziemy musieli dokonać
| |
| analogicznej permutacji wektora prawej strony
| |
|
| |
| Wynikiem takiego algorytmu jest więc rozkład
| |
|
| |
| <center><math>\displaystyle PA = LU,
| |
| </math></center>
| |
|
| |
| gdzie <math>\displaystyle P</math> jest pewną (zerojedynkową) macierzą permutacji (tzn. macierzą
| |
| identyczności z przepermutowanymi wierszami).
| |
|
| |
| Oprócz wyboru elementu głównego w kolumnie, stosuje się czasem inne strategie,
| |
| m.in. wybór w kolumnie oraz tzw. ''wybór pełny'', gdy elementu głównego szukamy w
| |
| ''całej'' podmacierzy <math>\displaystyle A(k:N,k:N)</math>, co znacznie zwiększa liczbę porównań
| |
| niezbędnych do wskazania elementu głównego, ale też trochę poprawia własności
| |
| numeryczne takiego algorytmu.
| |
|
| |
| W praktyce całą informację o wykonanych permutacjach przechowujemy nie w
| |
| zerojedynkowej macierzy j.w. ale w jednym wektorze.
| |
|
| |
| {{algorytm|Rozkład LU z wyborem elementu głównego w kolumnie||
| |
| <pre>
| |
|
| |
| P <nowiki>=</nowiki> 1:N; /* tu zapamiętujemy wykonane permutacje */
| |
| for k<nowiki>=</nowiki>1:N-1
| |
|
| |
| w wektorze A(k:N,k) znajdź element główny <math>\displaystyle a_{pk}</math>;
| |
| zamień ze sobą wiersze A(k,k:N) i A(p,k:N);
| |
| P(k) <nowiki>=</nowiki> p; P(p) <nowiki>=</nowiki> k;
| |
| if <math>\displaystyle a_{kk}</math>
| |
| STOP: macierz osobliwa!
| |
| end
| |
|
| |
| /* kontunuuj tak jak w algorytmie bez wyboru */
| |
|
| |
| for i<nowiki>=</nowiki>k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| |
| <math>\displaystyle a_{ik}</math> <nowiki>=</nowiki> <math>\displaystyle a_{ik}/a_{ii}</math>;
| |
| end
| |
| for j<nowiki>=</nowiki>k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
| |
| for i<nowiki>=</nowiki>k+1:N
| |
| <math>\displaystyle a_{ij}</math> -<nowiki>=</nowiki> <math>\displaystyle a_{ik}a_{kj}</math>;
| |
| end
| |
| end
| |
| end
| |
| </pre>}}
| |
|
| |
| Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
| |
| już rozkład LU wykonany z wyborem w kolumnie.
| |
|
| |
| {{algorytm|Rozwiązywanie układu równań metodą eliminacji Gaussa z wyborem elementu głównego w kolumnie||
| |
| <pre>
| |
|
| |
| znajdź rozkład <math>\displaystyle PA = LU</math>;
| |
| rozwiąż względem <math>\displaystyle y</math> układ z macierzą górną trójkątną <math>\displaystyle Uy = Pb</math>;
| |
| rozwiąż względem <math>\displaystyle x</math> układ z macierzą dolną trójkątną <math>\displaystyle Lx = y</math>;
| |
| </pre>}}
| |
|
| |
| <div class<nowiki>=</nowiki>"thumb tright"><div><flash>file<nowiki>=</nowiki>eliminacjagaussa.swf</flash><div.thumbcaption>Eliminacja Gaussa z wyborem elementu głównego</div></div></div>
| |
|
| |
| Dla niektórych ważnych klas macierzy wiadomo, że rozkład LU jest wykonalny ''bez wyznaczania elementu głównego'', co istotnie może zmniejszyć całkowity czas
| |
| działania algorytmu. Jest tak m.in. dla macierzy
| |
| * symetrycznych, dodatnio określonych: <math>\displaystyle A=A^T</math> oraz <math>\displaystyle x^TAx > 0</math>, <math>\displaystyle \forall x
| |
| \neq 0</math>,
| |
| * silnie diagonalnie dominujących: macierz <math>\displaystyle A</math> (lub <math>\displaystyle A^T</math>) spełnia
| |
|
| |
| <center><math>\displaystyle |a_{ii}| > \sum_{j\neq i} |a_{ij}|, \qquad \forall i.
| |
| </math></center>
| |
|
| |
| ==Pamięć hierarchiczna komputerów. Biblioteki BLAS, LAPACK==
| |
|
| |
| W poprzednim wykładzie sprawdziliśmy, jaki jest koszt obliczeniowy algorytmu
| |
| eliminacji Gaussa. Jednak w obecnym świecie istotna jest nie tylko liczba
| |
| operacji arytmetycznych, ale także koszt pobrania danych z pamięci. Za chwilę
| |
| zobaczymy, że poprzez ''reorganizację kolejności obliczeń'' w algorytmie
| |
| eliminacji Gaussa, możemy dostać algorytm (tzw. algorytm blokowy), którego
| |
| implementacja, choć znacznie mniej czytelna niż powyżej, będzie ''znacznie''
| |
| szybsza!
| |
|
| |
| Bez dostatecznie szybkiej pamięci, procesor -- zamiast liczyć -- będzie
| |
| większość czasu czekał na dane, a jego efektywność drastycznie spadnie. Z
| |
| niewielką przesadą można powiedzieć, że
| |
|
| |
| <blockquote> W optymalizacji szybkości działania programu numerycznego,
| |
| obecnie cała walka idzie o to, by procesor przez cały czas ''miał co liczyć''.
| |
| </blockquote>
| |
|
| |
| Szczególnie jest to widoczne w algorytmach, które wykonują bardzo dużo operacji
| |
| na dużej liczbie
| |
| danych --- a tak jest m.in. w algorytmach algebry liniowej takich jak mnożenie
| |
| dwóch macierzy, czy rozwiązywanie układów równań liniowych: te algorytmy
| |
| najczęściej operują na <math>\displaystyle O(N^2)</math> danych i wykonują aż <math>\displaystyle O(N^3)</math> działań.
| |
|
| |
| Ponieważ szybka pamięć komputerowa jest jednocześnie bardzo droga,
| |
| konstruktorzy komputerów osobistych stoją przed wykluczającymi się celami: z
| |
| jednej strony powinni zapewnić użytkownikowi jak najwięcej pamięci, z drugiej
| |
| zaś -- chcieliby zapewnić użytkownikowi jak najszybszą pamięć. Sprzeczność tych
| |
| dążeń ujawnia się, gdy dołączyć do nich trzecie, ale zasadnicze wymaganie:
| |
| całość ma być w rozsądnej cenie... Ostatecznie więc, z biegiem lat dramatycznie
| |
| pogłębia się przepaść pomiędzy prędkością (podwajającą się, zgodnie z
| |
| heurystycznym prawem Moore'a, co półtora roku) procesora, a prędkością pamięci
| |
| RAM, do której procesor musi się odwoływać.
| |
|
| |
| Powszechnie przyjętym sposobem pogodzenia tych sprzeczności jest pamięć
| |
| hierarchiczna. Koncept polega na tym, że część pamięci, która najczęściej komunikuje się z procesorem,
| |
| jest bardzo szybka (lecz jest jej relatywnie mało), natomiast pozostała pamięć
| |
| jest wolniejsza, za to może jej być bardzo dużo.
| |
|
| |
| W praktycznej realizacji, zamiast dwóch poziomów pamięci (mała-szybka i
| |
| duża-wolna), w komputerze występuje wiele poziomów:
| |
| * rejestry procesora
| |
| * ''cache ''(pamięć podręczna) procesora
| |
| * ''cache ''drugiego poziomu (ostatnio także wbudowywana do procesora)
| |
| * pamięć operacyjna (RAM): główna pamięć komputera
| |
| * pamięć zewnętrzna (np. twardy dysk)
| |
| * pamięć masowa (CD-ROM, taśma streamera, itp.)
| |
|
| |
| Efektywność działania komputera, zwłaszcza w wielkoskalowych obliczeniach
| |
| numerycznych, bardzo istotnie zależy od skutecznego wykorzystania hierarchii
| |
| pamięci tak, by jak największa część operacji była wykonywana na zmiennych
| |
| znajdujących się w danej chwili w jak najszybszej pamięci procesora.
| |
|
| |
| Kluczem do skutecznego wykorzystania hierarchii pamięci jest zasada
| |
| lokalności w czasie i w przestrzeni:
| |
|
| |
| <blockquote>
| |
| * '''Lokalność w czasie:''' Używać danego fragmentu pamięci intensywnie, ale rzadko.
| |
| * '''Lokalność w przestrzeni (adresowej):''' W danej chwili, odnosić
| |
| się do adresów pamięci leżących blisko siebie.
| |
|
| |
| </blockquote>
| |
|
| |
| Zachowanie zasady lokalności w czasie jest bardzo ważne dla efektywnego
| |
| wykorzystania cache'a -- ze względu na małe rozmiary tej pamięci, wymiana
| |
| zmiennych może być tam intensywna. Zasada lokalności w przestrzeni też jest
| |
| ważna, przy czym nie tylko ze względu na efektywne wykorzystanie cache'a, ale
| |
| także dla efektywnego wykorzystania pamięci
| |
| wirtualnej.
| |
|
| |
| ===Jak napisać kod źle wykorzystujący pamięć podręczną?===
| |
|
| |
| Choć programista nie ma bezpośredniego wpływu na opisane powyżej działania
| |
| systemu operacyjnego i ''hardware '''u (zarządzających, odpowiednio, pamięcią
| |
| wirtualną i ''cache ''), to przez właściwe projektowanie algorytmów -- a
| |
| zwłaszcza: ich ''właściwą'' implementację -- może spowodować, że jego
| |
| programy nie będą zmuszały komputera do nieracjonalnych zachowań. Jak się
| |
| okazuje, całkiem łatwo nieświadomie tworzyć programy przeczące zasadom
| |
| efektywnego wykorzystania hierarchii pamięci. Na potwierdzenie zacytujmy za
| |
| Dongarrą \cite{Dongarra} klasyczny już przykład mnożenia dwóch macierzy.
| |
|
| |
| W programie wykonujemy operację mnożenia dwóch macierzy <math>\displaystyle 1024\times 1024</math> przy
| |
| użyciu kilku ''matematycznie równoważnych'' algorytmów (nazwaliśmy je umownie
| |
| ijk, ikj, bikj(<math>\displaystyle \cdot</math>) --- nazwy pochodzą od sposobu organizacji pętli, zob.
| |
| poniżej), zaimplementowanych w programie C, wykorzystującym technikę
| |
| pozwalającą przechowywać macierze w pamięci komputera kolumnowo (zob.
| |
| Rozdział [[##sec:macierze-w-komputerze|Uzupelnic: sec:macierze-w-komputerze ]]). Dla porównania zmierzyliśmy czas
| |
| wykonania tej samej operacji przy użyciu wyspecjalizowanych bibliotek z
| |
| pakietów BLAS (algorytm DGEMM) i ATLAS (algorytm ATLAS DGEMM) --- zob. także
| |
| Rozdział [[##sec:blaslapack|Uzupelnic: sec:blaslapack ]]. Oto jakie wyniki uzyskaliśmy dla obliczeń w
| |
| arytmetyce podwójnej precyzji <code>double</code> na maszynie z procesorem AMD Duron
| |
| i zegarem 1.1 GHz:
| |
|
| |
| {| border=1
| |
| |+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
| |
| |-
| |
| |
| |
| Algorytm || ijk || ikj || bikj(16) || bikj(32) || DGEMM || ATLAS DGEMM
| |
| |-
| |
| |
| |
| Czas (s) || 320.49 || 24.28 || 8.68 || 30.45 || 25.72 || 2.58
| |
| |-
| |
| |
| |
| Mflop/s || 10.06 || 132.67 || 371.11 || 105.79 || 125.24 || 1248.53
| |
| |-
| |
| |
| |
|
| |
| |}
| |
|
| |
| Jak widać, różnice pomiędzy --- podkreślmy, matematycznie równoważnymi ---
| |
| algorytmami są bardzo znaczące; w szczególności algorytm ijk wydaje się nie do
| |
| przyjęcia! Powodem różnic musi być, ponieważ liczba wykonanych operacji
| |
| arytmetycznych jest identyczna, odmienne wykorzystanie pamięci ''cache ''
| |
| wynikające z organizacji dostępu do pamięci w naszych algorytmach.
| |
| Przedyskutujmy to dokładniej.
| |
|
| |
| ====Algorytm ijk====
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
| [ijk]
| |
| /* ijk */
| |
| for (i <nowiki>=</nowiki> 0; i < N; i++)
| |
| for (j <nowiki>=</nowiki> 0; j < N; j++)
| |
| for (k <nowiki>=</nowiki> 0; k < N; k++)
| |
| C[i*N+j] +<nowiki>=</nowiki> A[i*N+k]*B[k*N+j];
| |
| </pre></div>
| |
|
| |
| Jest to algorytm, który zapewne większości z Czytelników pierwszy przyszedłby
| |
| do głowy, gdyż realizuje wprost powszechnie znaną regułę mnożenia macierzy
| |
| "wiersz przez kolumnę". W pamięci ''cache ''L1 mieści się 64KB danych i
| |
| jest ona podzielona na 512 klas odpowiadających kolejnym liniom pamięci. W
| |
| każdej klasie można zmieścić 2 linie pamięci (Duron ma
| |
| ''2-way set associative cache''), a w każdej linia pamięci (i ''cache '''a) składa się z 64
| |
| bajtów, czyli mieści 8 liczb <code>double</code>.
| |
|
| |
| Odwołując się w najgłębszej pętli do kolejnych elementów macierzy <math>\displaystyle A</math> ''oraz'' <math>\displaystyle B</math> powodujemy, że przy odwoływaniu się do <math>\displaystyle B</math>,
| |
| ''cache miss ''następuje praktycznie w każdym kroku. Dzieje się tak dlatego,
| |
| że wymiary naszych macierzy są wielokrotnością 512: odwołując się
| |
| do kolejnych <code>B[k*N+j]</code>, <code>k</code> <nowiki>=</nowiki> <math>\displaystyle 0\ldots N</math>, odwołujemy się do co
| |
| 1024. elementu wektora B, a zatem kolejne odwołania lądują w tej samej sekcji
| |
| ''cache '''a -- która mieści ledwie 2 linie. Skutkiem tego, że w każdym
| |
| obrocie pętli mamy chybienie, jest złudzenie pracowania z komputerem ''bez'' pamięci ''cache ''(a nawet gorzej, bo ''cache miss ''dodatkowo
| |
| kosztuje) czyli jakby z komputerem o szybkości 10 MHz <nowiki>=</nowiki> 100 MHz/10 (bo
| |
| magistrala (''bus '') jest taktowana 100 MHz, a odwołanie do pamięci RAM
| |
| kosztuje mniej więcej 10 cykli magistrali). I rzeczywiście, wyniki zdają się to
| |
| potwierdzać.
| |
|
| |
| ====Algorytm ikj====
| |
|
| |
| Różni się on od poprzedniego jedynie
| |
| kolejnością dwóch wewnętrznych pętli:
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
| [ikj]
| |
| /* ikj */
| |
| for (i <nowiki>=</nowiki> 0; i < N; i++)
| |
| for (k <nowiki>=</nowiki> 0; k < N; k++)
| |
| for (j <nowiki>=</nowiki> 0; j < N; j++)
| |
| C[i*N+j] +<nowiki>=</nowiki> A[i*N+k]*B[k*N+j];
| |
| </pre></div>
| |
|
| |
| Okazuje się, że taka prosta zmiana dramatycznie poprawia sytuację!
| |
|
| |
| Tym razem, w odwołaniu do <math>\displaystyle B</math> w wewnętrznej pętli, ''cache miss '' będzie
| |
| następował ośmiokrotnie rzadziej, gdyż odwołując się do ''kolejnych''
| |
| elementów wektora <code>B</code>, znacznie częściej odwołujemy się do danych
| |
| znajdujących się w ''cache '',
| |
| zachowując zasadę ''lokalności w przestrzeni'': ponieważ w linii ''cache '''a
| |
| mieści się osiem ''kolejnych'' elementów wektora <code>B</code>. Stąd
| |
| znaczące przyspieszenie (więcej niż ośmiokrotne --- dlaczego?).
| |
|
| |
| ====Algorytm bikj()====
| |
|
| |
| Algorytm bikj(16) jest prostą wersją algorytmu blokowego operującego w sposób
| |
| "ikj" na blokach macierzy wymiaru <math>\displaystyle 16\times 16</math>:
| |
|
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
| [bikj(16)]
| |
| /* bikj(16) */
| |
| for (i <nowiki>=</nowiki> 0; i < N; i+<nowiki>=</nowiki>16)
| |
| for (k <nowiki>=</nowiki> 0; k < N; k+<nowiki>=</nowiki>16)
| |
| for (j <nowiki>=</nowiki> 0; j < N; j+<nowiki>=</nowiki>16)
| |
| for (ii <nowiki>=</nowiki> i; ii < i+15; ii++)
| |
| for (kk <nowiki>=</nowiki> k; kk < k+15; kk++)
| |
| for (jj <nowiki>=</nowiki> j; jj < j+15; jj++)
| |
| C[ii*N+jj] +<nowiki>=</nowiki> A[ii*N+kk]*B[kk*N+jj];
| |
|
| |
| </pre></div>
| |
|
| |
| (podobnie działał algorytm bikj(32), tym razem na blokach <math>\displaystyle 32\times 32</math>).
| |
|
| |
| Wadą poprzedniego algorytmu, ikj, było to, że w dwu zewnętrznych pętlach powracał do
| |
| wszystkich <math>\displaystyle N^2</math> wartości <math>\displaystyle C</math> i <math>\displaystyle B</math>, przecząc zasadzie ''lokalności w
| |
| czasie''. Wydzielenie w algorytmie bijk(16) operacji na blokach macierzy ma na celu uniknięcie tego
| |
| problemu: trzy najbardziej zewnętrzne pętle to <math>\displaystyle (1024/16)^3</math> obrotów, czyli
| |
| czterokrotnie mniej niż dwie najbardziej zewnętrzne pętle w poprzednim
| |
| algorytmie. Gdyby udało się zachować liczbę ''cache misses ''na poprzednim poziomie,
| |
| można byłoby liczyć na czterokrotne przyspieszenie. ('''Do sprawdzenia: I tak z grubsza jest:
| |
| teoretycznie wszystkie 3 podmacierze powinny mieścić się w cache'u.''')
| |
|
| |
| ====Algorytmy DGEMM i ATLAS DGEMM====
| |
|
| |
| Algorytm DGEMM to był algorytm mnożenia macierzy z pakietu BLAS --- jest to
| |
| właśnie profesjonalny algorytm blokowy, ale niezoptymalizowany na naszą
| |
| architekturę. Dopiero algorytm DGEMM podrasowany w pakiecie ATLAS dał nam
| |
| sprawność wynoszącą 1.2 Gflopów na sekundę -- czyli praktycznie
| |
| maksimum (''teoretycznie'', z Durona można wycisnąć dwa flopy w cyklu
| |
| zegara, co dawałoby <math>\displaystyle r_{\max}</math> <nowiki>=</nowiki> 2.2 Gflop/s, ale w praktyce jest to mało
| |
| prawdopodobne) tego,
| |
| co da się wycisnąć z tej wersji Durona na liczbach podwójnej precyzji.
| |
|
| |
| ===Macierze w pamięci komputera===
| |
|
| |
| Do implementacji algorytmów numerycznych powszechnie używa się dwóch języków:
| |
| Fortranu i C. Zgodnie z naszą konwencją, skupimy się na programach w C, nie
| |
| możemy jednak przejść obojętnie wobec faktu, że jest bardzo wiele znakomitego
| |
| oprogramowania numerycznego w Fortranie. W Rozdziale [[##sec:FortranC|Uzupelnic: sec:FortranC ]] zajmiemy się
| |
| metodą włączenia podprogramów fortranowskich do programu w C. Zanim jednak to
| |
| uczynimy, musimy zauważyć, że oba języki przechowują macierze w pamięci
| |
| komputera w odmienny sposób, co stanowi źródło potencjalnych poważnych kłopotów
| |
| na styku tych języków. Dlatego w niniejszym rozdziale zechcemy szczegółowo
| |
| przyjrzeć się temu, jak oba te języki przechowują w pamięci macierze i
| |
| opracować sposoby postępowania gwarantujące zgodność programów w obu
| |
| językach.
| |
|
| |
| W Fortranie, elementy macierzy są przechowywane w pamięci kolumnami, tzn. jeśli
| |
| mamy do czynienia z macierzą prostokątną <math>\displaystyle n\times m</math> o elementach <math>\displaystyle a_{ij}</math>,
| |
| <math>\displaystyle i=1\ldots n</math>, <math>\displaystyle j=1\ldots m</math>,
| |
|
| |
| <center><math>\displaystyle
| |
| \begin{pmatrix}
| |
| a_{11} & \cdots & a_{1m}\\
| |
| \vdots & & \vdots\\
| |
| a_{n1} & \cdots & a_{nm}
| |
| \end{pmatrix} .
| |
| </math></center>
| |
|
| |
| to kolejne miejsca w przestrzeni adresowej
| |
| zajmują elementy
| |
| <center><math>\displaystyle a_{11}, a_{21},\ldots, a_{n1}, a_{12}, a_{22}, \ldots, a_{n2},
| |
| \ldots a_{nm}.
| |
| </math></center>
| |
|
| |
| Dla odmiany, C przechowuje w pamięci elementy macierzy wierszami, tzn. kolejno
| |
| <center><math>\displaystyle a_{11}, a_{12},\ldots, a_{1m}, a_{21}, a_{22}, \ldots, a_{2m}, \ldots
| |
| a_{nm}.</math></center>
| |
|
| |
| Co więcej, standard nie gwarantuje, że kolejne wiersze macierzy będą
| |
| przechowywane w przylegających do siebie obszarach pamięci. Bierze się to stąd,
| |
| że w C macierz dwuwymiarowa jest w istocie tablicą wskaźników do kolejnych
| |
| wierszy.
| |
|
| |
| To zaś powoduje kolejne komplikacje. Otóż w procedurach numerycznych bardzo
| |
| często pragniemy dynamicznie przydzielać pamięć na tablice, to znaczy dopiero
| |
| w trakcie działania procedury, gdyż dopiero w trakcie obliczeń procedura
| |
| otrzymuje np. informację o rozmiarze potrzebnej macierzy. Przykładowo, program w C,
| |
| który wykonywałby dynamicznie przydzielałby pamięć na dwuwymiarową tablicę,
| |
| musiałby:
| |
| * przydzielić pamięć na tablicę wskaźników do wierszy
| |
| * każdemu wskaźnikowi do wiersza przydzielić pamięć na pojedynczy wiersz
| |
|
| |
| To jest jeden z licznych powodów, dla których posługując się macierzami w C
| |
| będziemy stosowali pewien prosty ''trick ''.
| |
|
| |
| Otóż przyjmiemy konwencję, że nie będziemy wcale korzystać z macierzy
| |
| dwu- i więcejwymiarowych, a elementy macierzy zapiszemy do jednego odpowiednio
| |
| długiego wektora. W ten sposób wszystkie elementy macierzy zajmą spójny
| |
| obszar pamięci. Przykładowo, macierz wymiaru <math>\displaystyle n\times m</math> będziemy zapisywali do wektora
| |
| o długości <math>\displaystyle n\cdot m</math>.
| |
|
| |
| Mając pełną dowolność wyboru sposobu ułożenia elementów macierzy w wektorze,
| |
| wybierzemy zazwyczaj układ kolumnowy (czyli fortranowski) (niektóre biblioteki w
| |
| C (np. [http://www.fftw.org FFTW]) wymagają jednak układu wierszowego!),
| |
| co ma dwie zalety: po pierwsze, da nam zgodność z Fortranem, a po drugie, może
| |
| potencjalnie uprościć nam optymalizację "książkowych" algorytmów, których
| |
| większość była i często nadal jest pisana z myślą o implementacji w Fortranie.
| |
| Złudzenie korzystania z macierzy dwuwymiarowej da nam zaś makro, lokalizujące
| |
| <math>\displaystyle (i,j)</math>-ty element macierzy w wektorze przechowującym jej elementy. Dodatkowo,
| |
| makro tak skonstruowaliśmy, aby móc indeksować elementy macierzy poczynając od
| |
| 1, czyli <math>\displaystyle a_{ij}</math>, <math>\displaystyle i=1\ldots n</math>, <math>\displaystyle j=1\ldots m</math>.
| |
|
| |
| Poniżej pokazuje przykładowy fragment kodu realizującego opisaną wyżej ideę.
| |
|
| |
| Zwróćmy uwagę na dwa sposoby odwoływania się do elementów macierzy. Za pierwszym
| |
| razem, odwołujemy się do kolejnych elementów wektora <code>matrix</code>, gdyż pętle są
| |
| ustawione tak, by przechodzić przez macierz wzdłuż kolejnych kolumn. Dlatego nie
| |
| jest tu konieczne użycie makra <code>IJ()</code>, a sprytne wykorzystanie
| |
| pointera <code>ptr</code>
| |
| pozwala na zrezygnowanie z operacji mnożenia przy odwoływaniu się do kolejnych
| |
| elementów macierzy.
| |
|
| |
| Za drugim razem chcemy wydrukować naszą macierz na ekranie, musimy więc
| |
| odwoływać się do kolejnych ''wierszy'' macierzy (a więc, z punktu
| |
| wykorzystania cache i hierarchii pamięci, fatalnie! -- na szczęście I/O jest
| |
| znacznie wolniejszy od najwolniejszej nawet pamięci RAM). Tym razem więc nie
| |
| unikniemy wywołania makra <code>IJ()</code> (i obliczania wyrażenia <code>i+j*N</code>) przy
| |
| każdym obrocie wewnętrznej pętli.
| |
|
| |
| Inne zalety takiego podejścia do przechowywania macierzy w C to:
| |
| * łatwy dostęp do takich macierzy z funkcji fortranowskich
| |
| * właściwie opracowane makro <code>IJ()</code> pozwala na ominięcie
| |
| problemu indeksowania elementów macierzy (w C macierze indeksujemy od zera, gdy
| |
| tymczasem we wszelkich "ludzkich" algorytmach (i, dodajmy, np. w Octave i
| |
| MATLABie), elementy macierzy indeksowane są od "1";
| |
| * jawny sposób ułożenia elementów macierzy w pamięci zwiększa przenośność i
| |
| odporność implementowanych procedur
| |
|
| |
| Ceną jaką za to płacimy, jest używanie za każdym razem makra w celu odwołania
| |
| się do konkretnego elementu macierzy. Ponadto, można byłoby się przyczepić do
| |
| niepotrzebnie wielokrotnego wyznaczania tego samego iloczynu <code>j*N</code>, gdy
| |
| odwołujemy się do kolejnych elementów kolumny macierzy. Jest to (niewygórowana,
| |
| moim zdaniem) cena, jaką płacimy za przejrzystość, czytelność, elastyczność i
| |
| "wielojęzyczność" (C/Fortran) naszego programu. Dodajmy, że opisane podejście
| |
| nie jest niczym nowym w praktyce numerycznej i jest stosowane w wielu dobrych
| |
| bibliotekach numerycznych; można tu wymienić np. doskonały skądinąd pakiet
| |
| CVODE (macierz w wektorze plus makra <code>IJ()</code>) czy też pakiet CLAPACK
| |
| (macierz w wektorze), zob.
| |
| \cite{clapack-howto}).
| |
|
| |
| Przypomnijmy przy okazji, że ze względu na konstrukcję ''cache '''a spotykaną
| |
| np. w procesorach Intela i AMD, czasem warto stosować tzw. ''array padding ''
| |
| w przypadku, gdy mamy do czynienia z macierzami o wymiarach będących dużą
| |
| potęgą dwójki, zob. Rozdział [[##sec:cache:example|Uzupelnic: sec:cache:example ]].
| |
