|
|
(Nie pokazano 77 wersji utworzonych przez 3 użytkowników) |
Linia 1: |
Linia 1: |
| ''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi� poprawki''
| |
|
| |
|
| ==Rozwiązywanie układów równań liniowych== | | <!-- |
| | Konwertowane z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/ pawlik1/latex2mediawiki.php. |
| | Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki |
| | wprowadzi� przykry@mimuw.edu.pl |
| | --> |
| | |
| | =Własności zadania obliczeniowego i algorytmu numerycznego= |
|
| |
|
| Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń
| | {{powrot |Metody numeryczne | do strony głównej |
| zmiennoprzecinkowych, przy czym wiele znanych, ''matematycznie równoważnych'' metod
| | przedmiotu <strong>Metody numeryczne</strong>}} |
| 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
| | ==Uwarunkowanie zadania obliczeniowego== |
| 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>.
| | Jak zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując |
| | zadanie obliczeniowe, są z natury rzeczy wartościami zaburzonymi. Źródłem tych |
| | zaburzeń są: |
| | * błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (np. 0.1 nie jest równe dokładnie <math>1/10</math>) |
| | * błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie <math>f(x) = a</math>, ale <math>a</math> jest rezultatem innej symulacji), a także |
| | * błędy pomiarowe w zadaniach praktycznych (np. chcemy policzyć numeryczną prognozę pogody, ale temperaturę wyjściową znamy tylko z dokładnością do 0.1 stopnia, co gorsza --- z niektórych stacji w ogóle nie mamy danych) |
| | |
| | Okazuje się, że powszechna intuicja, że małe zaburzenia danych powinny dawać |
| | małe zaburzenia wyniku, nie znajduje potwierdzenia nawet w bardzo prostych |
| | przypadkach. Z drugiej strony, umiejętność oceny jakościowego <strong>wpływu |
| | zaburzenia danych na wynik</strong> jest kapitalna w świecie obliczeń numerycznych w |
| | ogólności, a w szególności --- inżynierskich. |
|
| |
|
| W
| | Wprowadza się pojęcie <strong>uwarunkowania</strong> zadania, to znaczy jego podatności na |
| praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
| | zaburzenia danych. Dla przejrzystości przypuśćmy, że nasze zadanie obliczeniowe |
| specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
| | polega na wyznaczeniu <math>f(x)</math> dla danego <math>x</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ń
| | Jak bardzo będzie odległe |
| liniowych, takie jak:
| | <math>f(\widetilde{x})</math>, gdy <math>\widetilde{x}\approx x</math>? Rozważa się dwa przypadki: |
| * metoda wyznacznikowa (wzory Cramera) | | * <strong>uwarunkowanie względne</strong>: jak względne zaburzenie danych wpływa na błąd względny wyniku: <center><math>\frac{||f(x) - f(\widetilde{x})||}{||f(x)||} \leq \mbox{cond} _{rel}(f,x) \cdot \frac{||x - \widetilde{x}||}{||x||}</math></center> |
| * obliczenie macierzy <math>\displaystyle A^{-1}</math> i następnie <math>\displaystyle x = A^{-1}b</math> | | Najmniejszy mnożnik <math>\mbox{cond} _{rel}(f,x)</math> spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>. |
| | * <strong>uwarunkowanie bezwzględne</strong>: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku: <center><math>||f(x) - f(\widetilde{x})|| \leq \mbox{cond} _{abs}(f,x) \cdot ||x - \widetilde{x}||</math></center> |
| | Najmniejszy mnożnik <math>\mbox{cond} _{abs}(f,x)</math> spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>. |
| | |
| | Powiemy, że zadanie <math>f(x)</math> jest |
| | * <strong>dobrze uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \approx 1</math>, |
| | * <strong>źle uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \gg 1</math>, |
| | * <strong>źle postawione</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) = +\infty</math>. |
| | |
| | Teraz już rozumiemy, że redukcja cyfr przy odejmowaniu jest tylko |
| | odzwierciedleniem faktu, że zadanie odejmowania dwóch bliskich liczb jest po |
| | prostu zadaniem źle uwarunkowanym! |
|
| |
|
| '''nie nadają się''' do numerycznego rozwiązywania takich zadań.
| | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
| | <span style="font-variant:small-caps;">Przykład: Uwarunkowanie zadania obliczenia sumy</span> |
| | <div class="solution" style="margin-left,margin-right:3em;"> |
|
| |
|
| O tym, jak ''skutecznie'' rozwiązywać takie zadania, jakie jest ich
| | Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia <math>s(x,y) = x + y</math> ma |
| uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z
| | <center><math> |
| dobrych metod jest metoda eliminacji Gaussa.
| | \mbox{cond} _{abs}(s, (a,b)) = 1, \qquad \mbox{cond} _{rel}(s, (a,b)) = \frac{|a|+|b|}{|a+b|} |
| | |
| ===Proste układy równań===
| |
| | |
| Niektóre układy równań można bardzo łatwo rozwiązać. Zgodnie z zasadą
| |
| | |
| Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
| |
| | |
| 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"?
| |
| | |
| ''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> | | </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>.
| | Tak więc, gdy <math>a\approx -b</math>, to <math>\mbox{cond} _{rel}(s, (a,b)) \approx +\infty</math> i zadanie |
| | | jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może |
| W
| | skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego, |
| praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
| | najczęściej rzeczywiście tak będzie... |
| specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
| | </div></div> |
| | |
| 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ą
| |
| | |
| Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
| |
| | |
| 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ą==== | | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
| | <span style="font-variant:small-caps;">Przykład</span> |
| | <div class="solution" style="margin-left,margin-right:3em;"> |
|
| |
|
| Rozważmy układ z macierzą
| | Dla różniczkowalnej funkcji skalarnej <math>f : R \rightarrow R</math> mamy |
| trójkątną <math>\displaystyle A</math>. Będą nas szczególnie interesować macierze
| | <center><math> |
| ''trójkątne górne'', dla których <math>\displaystyle a_{i,j}=0</math> gdy <math>\displaystyle i>j</math>, oraz
| | |f(x) - f(\widetilde{x})| \approx |f'(x) | | x - \widetilde{x} | |
| 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> | | </math></center> |
|
| |
|
| Układ z nieosobliwą macierzą trójkątną górną
| | i w konsekwencji dla zadania obliczenia <math>f(x)</math> dla danego <math>x</math> mamy, przy |
| | założeniu małych zaburzeń, |
| | <center><math> |
| | \mbox{cond} _{abs}( f, x) = |f'(x)|, \qquad \mbox{cond} _{rel}( f, x) = |
| | \frac{|f'(x)|\cdot|x|}{|f(x)|}</math></center> |
|
| |
|
| <center><math>\displaystyle
| | </div></div> |
| U\,\vec x\;=\;\vec c,
| |
| </math></center> | |
|
| |
|
| <math>\displaystyle U=(u_{i,j})</math>, <math>\displaystyle \vec c=(c_j)</math>, można rozwiązać stosując algorytm:
| | Jest zrozumiałe, że złe uwarunkowanie jest szkodliwe w praktyce numerycznej: zaburzenia danych są nieuniknione, a źle uwarunkowane zadanie tylko je wzmocni na wyjściu. Jednak za jakiś czas zobaczysz wyjątkową sytuację, gdy [[MN13#Odwrotna metoda potęgowa|złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, ale wręcz <strong>pomaga</strong>]] szybciej rozwiązać zadanie główne! |
|
| |
|
| {{algorytm||
| | ==Rozkład algorytmu względem informacji== |
| <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>;
| |
| }}
| |
|
| |
|
| (Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy
| | <strong>Algorytm</strong> to dokładnie określona i dozwolona w danym modelu |
| implikuje, że <math>\displaystyle u_{i,i}\ne 0</math>, <math>\displaystyle \forall i</math>.) Podobnie, układ
| | obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego |
| <math>\displaystyle L\vec x=\vec c</math> rozwiązujemy algorytmem:
| | zadania (w sposób dokładny lub przybliżony). |
|
| |
|
| {{algorytm||
| | Z każdym algorytmem związany jest operator |
| <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>;
| |
| }}
| |
|
| |
|
| Oba algorytmy wymagają rzędu <math>\displaystyle N^2/2</math> mnożeń lub dzieleń i
| | <center><math>{\bf ALG}:\,F\longrightarrow G, |
| <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\vec x = \vec b,
| |
| </math></center> | | </math></center> |
|
| |
|
| gdy <math>\displaystyle Q</math> jest macierzą ortogonalną, to znaczy <math>\displaystyle Q^TQ = I</math>. Rzeczywiście, z
| | taki że <math>{\bf ALG}(f)</math> jest wynikiem działania algorytmu |
| ortogonalności mamy natychmiast, że
| | w arytmetyce idealnej dla danej <math>f</math>. |
|
| |
|
| <center><math>\displaystyle | | Zauważmy, że wynik <math>{\bf ALG}(f)</math> działania algorytmu nie |
| \vec x = Q^T \vec b | | zależy bezpośrednio od <math>f</math>, ale raczej od <strong>informacji</strong> |
| </math></center> | | o <math>f</math> (uzyskanej dzięki poleceniu <math>{\cal IN}</math>). Informacja |
| | ta może być <strong>pełna</strong> albo tylko <strong>częściowa</strong>. |
| | Informacja jest pełna gdy, np. |
| | <math>f=(f_1,\ldots,f_n)\in R^n</math> i wczytamy wszystkie |
| | współrzędne <math>f_i</math>. Informacja może być częściowa, gdy |
| | <math>f</math> jest funkcją. Wtedy wiele danych może posiadać tę |
| | samą informację, co łatwo zaobserwować na przykładzie |
| | zadania całkowania. |
|
| |
|
| i w konsekwencji <math>\displaystyle x</math> można wyznaczyć kosztem takim jak koszt mnożenia macierzy
| | Niech <math>N:F\to \cup_{n=0}^\infty R^n</math> będzie |
| przez wektor, czyli <math>\displaystyle O(N^2)</math> operacji.
| | <strong>operatorem informacji</strong>, tzn. |
|
| |
|
| Podobnie, gdy <math>\displaystyle Q\in C^{N\times N}</math> jest unitarna, to znaczy <math>\displaystyle Q^HQ = I</math>,
| | <center><math>N(f)\,=\,(y_1,y_2,\ldots,y_n) |
| rozwiązaniem układu równań jest
| |
| | |
| <center><math>\displaystyle
| |
| \vec x = Q^H \vec b.
| |
| </math></center> | | </math></center> |
|
| |
|
| ===Metoda eliminacji Gaussa===
| | jest informacją o <math>f</math> zebraną przy idealnej realizacji |
| | | algorytmu. Zauważmy, że informacja jest pełna, gdy <math>N</math> jest |
| W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym | | przekształceniem różnowartościowym, tzn. jeśli |
| rozwiązywania układu równań
| | <math>f_1\ne f_2</math> implikuje <math>N(f_1)\ne N(f_2)</math>. |
| | W przeciwnym przypadku mamy do czynienia z informacją |
| | częściową. |
|
| |
|
| <center><math>\displaystyle Ax=b</math></center> | | Każdy algorytm <math>{\bf ALG}</math> może być przedstawiony jako złożenie |
| | operatora informacji i pewnego operatora |
| | <math>\varphi:N(F)\to G</math> zdefiniowanego równością |
|
| |
|
| okazuje się popularna
| | <center><math>\varphi\left(N(f)\right)\,=\,{\bf ALG}(f). |
| 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> | | </math></center> |
|
| |
|
| a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:
| | Zauważmy, że w przypadku informacji częściowej zwykle nie |
| | | istnieje algorytm dający dokładne rozwiązanie zadania dla |
| {{algorytm|title<nowiki>=</nowiki>Algorytm eliminacji Gaussa|
| | każdej danej <math>f\in F</math>, ponieważ dla danych o tej samej |
| | | informacji mogą istnieć różne rozwiązania. |
| Znajdź rozkład <math>\displaystyle A=LU</math>;
| |
| Rozwiąż <math>\displaystyle Ly = b</math>;
| |
| Rozwiąż <math>\displaystyle Ux = y</math>;
| |
| }}
| |
| | |
| 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
| | ==Problem wyboru algorytmu== |
| \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
| | Wybór algorytmu jest najistotniejszą częścią całego procesu |
| * <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest | | numerycznego rozwiązywania zadania. Kierujemy się przy tym przede |
| kopią pierwszego wiersza <math>\displaystyle A</math>,
| | wszystkim następującymi kryteriami: |
| * <math>\displaystyle l_{21} = a_{21}/u_{11}</math>, więc pierwsza kolumna <math>\displaystyle L</math> powstaje przez | | * dokładnością algorytmu, |
| podzielenie wszytkich elementów wektora <math>\displaystyle a_{21}</math> przez element na diagonali
| | * złożonością algorytmu, |
| <math>\displaystyle a_{11}</math>, | | * własnościami numerycznymi algorytmu. |
| * <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
| | Przez dokładność algorytmu rozumiemy różnicę między |
| bloku <math>\displaystyle A_{22}</math> macierzy <math>\displaystyle A</math>,
| | rozwiązaniem dokładnym <math>S(f)</math> a rozwiązaniem |
| wymiaru <math>\displaystyle (N-1)\times (N-1)</math>.
| | <math>{\bf ALG}(f)</math> dawanym przez algorytm w arytmetyce idealnej. |
| | Jeśli <math>{\bf ALG}(f) = S(f)</math>, |
| | <math>\forall f \in F</math>, |
| | to algorytm nazywamy <strong>dokładnym</strong>. |
|
| |
|
| Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie
| | Mówiąc o złożoności, mamy na myśli złożoność pamięciową |
| rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
| | (zwykle jest to liczba stałych i zmiennych używanych przez |
| klasycznych pętli. Jest to ważne w praktyce numerycznej, gdyż rekurencja
| | algorytm), jak również złożoność obliczeniową. |
| kosztuje: zarówno pamięć, jak i czas.
| | Na złożoność obliczeniową algorytmu dla danej <math>f</math> składa |
| | się koszt uzyskania infomacji <math>y=N(f)</math> (zwykle jest on |
| | proporcjonalny do liczby wywołań polecenia <math>{\cal IN}</math>), oraz |
| | koszt <strong>kombinatoryczny</strong> przetworzenia tej informacji, aż do |
| | uzyskania wyniku <math>\varphi(y)</math>. Koszt kombinatoryczny zwykle |
| | mierzymy liczbą operacji arytmetycznych wykonywanych przez |
| | algorytm. |
|
| |
|
| Ponadto zauważmy, że opisany algorytm możemy wykonać w miejscu, nadpisując
| | Przez własności numeryczne algorytmu rozumiemy jego |
| 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
| | własności przy realizacji w arytmetyce <math>fl_\nu</math>. Temu |
| pamiętać, bo wiemy a priori, że tam są).
| | ważnemu tematowi poświęcimy teraz osobny paragraf. |
|
| |
|
| {{algorytm|title<nowiki>=</nowiki>Algorytm rozkładu LU|
| | ==Numeryczna poprawność algorytmu== |
|
| |
|
| for k<nowiki>=</nowiki>1:N-1
| | Pożądane jest, aby algorytm dawał "dobry" wynik zarówno |
| if <math>\displaystyle a_{kk}</math> <nowiki>=</nowiki><nowiki>=</nowiki> 0
| | w arytmetyce idealnej, jak i w arytmetyce <math>fl_\nu</math>. Niestety, |
| STOP;
| | jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm |
| end
| | jest dokładny, to w wyniku jego realizacji w <math>fl_\nu</math> możemy |
| for i<nowiki>=</nowiki>k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| | otrzymać wynik <math>fl_\nu({\bf ALG}(f))</math> daleko odbiegający od |
| <math>\displaystyle a_{ik}</math> <nowiki>=</nowiki> <math>\displaystyle a_{ik}/a_{ii}</math>; | | <math>S(f)</math>. W szczególności, prawie zawsze mamy |
| 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
| |
| }}
| |
|
| |
|
| Łatwo przekonać się, że <math>\displaystyle k</math>-ty obrót zewnętrznej pętli (tzn. <math>\displaystyle k</math>-ty krok
| | <center><math>S(f)\,\ne\,fl_\nu\left({\bf ALG}(f)\right). |
| 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ć {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> | | </math></center> |
|
| |
|
| jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od
| | Zauważmy również, że o ile z reguły znamy dokładne zachowanie |
| razu zetknie się z dzieleniem przez <math>\displaystyle a_{11}=0</math>. Ale wystarczy zamienić ze sobą
| | się algorytmu w arytmetyce idealnej dla danej informacji, to nie |
| kolejnością wiersze macierzy <math>\displaystyle A</math> (to znaczy, w układzie równań, zamienić ze sobą
| | można tego samego powiedzieć o jego zachowaniu się w arytmetyce |
| miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez
| | <math>fl_\nu</math>. W związku z tym powstaje pytanie, jak kontrolować błąd |
| problemu.
| | algorytmów wynikający z błędów zaokrągleń i jakie algorytmy |
| | uznamy za te o najwyższej jakości numerycznej. |
|
| |
|
| W praktyce obliczeniowej, aby uzyskać algorytm o {}{możliwie dobrych własnościach
| | Istnienie błędów reprezentacji liczb rzeczywistych powoduje, |
| numerycznych}, wykorzystujemy tzw. strategię ''wyboru elementu głównego w
| | że informacja <math>y=N(f)</math> o danej <math>f</math> nie jest w |
| kolumnie''.
| | ogólności reprezentowana dokładnie. Znaczy to, że zamiast na |
| Polega to na tym, że zanim wykonamy <math>\displaystyle k</math>-ty krok algorytmu rozkładu LU,
| | informacji dokładnej, dowolny algorytm będzie operować na |
| # szukamy w pierwszej kolumnie podmacierzy <math>\displaystyle A(k:N,k:N)</math> szukamy elementu o
| | informacji <strong>nieco zaburzonej</strong> <math>y_\nu</math>, tzn. zaburzonej na |
| największym module (taki element, na mocy założenia nieosobliwości macierzy,
| | poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm |
| jest niezerowy) --- to jest właśnie element główny
| | będzie w ogólności zaburzony na poziomie błędu reprezentacji. |
| # zamieniamy ze sobą wiersz pierwszy <math>\displaystyle A(k:N,k:N)</math> z wierszem, w którym
| | W najlepszym więc wypadku wynikiem działania algorytmu w <math>fl_\nu</math> |
| znajduje się element główny
| | będzie <math>(\varphi(y_\nu))_\nu</math> zamiast <math>\varphi(y)</math>. Algorytmy |
| # zapamiętujemy dokonaną permutację, bo potem --- gdy już przyjdzie do
| | dające tego rodzaju wyniki uznamy za posiadające najlepsze |
| rozwiązywania układu równań --- będziemy musieli dokonać
| | własności numeryczne w arytmetyce <math>fl_\nu</math> i nazwiemy numerycznie |
| analogicznej permutacji wektora prawej strony
| | poprawnymi. |
|
| |
|
| Wynikiem takiego algorytmu jest więc rozkład
| | Powiemy, że ciąg rzeczywisty |
| | <math>a_\nu=(a_{\nu,1},\ldots,a_{\nu,n})</math> |
| | (a właściwie rodzina ciągów <math>\{a_\nu\}_\nu</math>) jest |
| | <strong>nieco zaburzonym</strong> ciągiem <math>a=(a_1,\ldots,a_n)</math>, jeśli |
| | istnieje stała <math>K</math> taka, że dla wszystkich dostatecznie |
| | małych <math>\nu</math> zachodzi |
|
| |
|
| <center><math>\displaystyle PA = LU, | | <center><math> |
| | |a_{\nu,j} - a_j|\,\le\,K\,\nu\,|a_j|,\qquad 1\le j\le n, |
| </math></center> | | </math></center> |
|
| |
|
| gdzie <math>\displaystyle P</math> jest pewną (zerojedynkową) macierzą permutacji (tzn. macierzą
| | albo ogólniej |
| 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|title<nowiki>=</nowiki>Algorytm rozkładu LU z wyborem elementu głównego w kolumnie|
| | <center><math> |
| | \|a_\nu - a\|\,\le\,K\,\nu\,\|a\|</math>,</center> |
|
| |
|
| P <nowiki>=</nowiki> 1:N; /* tu zapamiętujemy wykonane permutacje */
| | gdzie <math>\|\cdot\|</math> jest pewną normą w <math>R^n</math>. W pierwszym |
| for k<nowiki>=</nowiki>1:N-1
| | przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim |
| | o zaburzeniu w normie <math>\|\cdot\|</math>. |
|
| |
|
| w wektorze A(k:N,k) znajdź element główny <math>\displaystyle a_{pk}</math>;
| | Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają |
| zamień ze sobą wiersze A(k,k:N) i A(p,k:N);
| | za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas |
| 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 */ | | <center><math>\|a_\nu - a\|_\infty \,=\, \max_{1\le j\le n} |a_{\nu,j} - a_j| |
| | \,\le\,K\,\nu\,\max_{1\le j\le n} |a_j|\,=\,K\,\nu\,\|a\|_\infty</math>,</center> |
|
| |
|
| for i<nowiki>=</nowiki>k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
| | i korzystając z faktu, że w przestrzeni skończenie wymiarowej |
| <math>\displaystyle a_{ik}</math> <nowiki>=</nowiki> <math>\displaystyle a_{ik}/a_{ii}</math>;
| | wszystkie normy są równoważne, otrzymujemy dla pewnych stałych |
| end
| | <math>K_1</math> i <math>K_2</math> |
| 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
| |
| }}
| |
|
| |
|
| Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
| | <center><math>\|a_\nu - a\|\,\le\,K_1\|a_\nu-a\|_\infty\,\le\, |
| już rozkład LU wykonany z wyborem w kolumnie.
| | K_1 K\,\nu\,\|a\|_\infty\,\le\,K_2 K_1 K\,\nu\,\|a\|</math>,</center> |
|
| |
|
| {{algorytm|title<nowiki>=</nowiki>Algorytm eliminacji Gaussa z wyborem elementu głównego w kolumnie|
| | czyli nierówność dla zaburzenia w normie, ze stałą <math>K = K_2 K_1 K</math>. |
|
| |
|
| znajdź rozkład <math>\displaystyle PA = LU</math>;
| | {{definicja|Algorytm numerycznie poprawny|Algorytm numerycznie poprawny| |
| 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>;
| |
| }}
| |
|
| |
|
| {}{Eliminacja Gaussa z wyborem elementu głównego} | | Algorytm <math>{\bf ALG}</math> rozwiązywania zadania |
| | nazywamy <strong>numerycznie poprawnym</strong> w zbiorze danych |
| | <math>F_0\subset F</math>, jeśli dla każdej danej <math>f\in F_0</math> |
| | wynik <math>fl_\nu({\bf ALG}(f))</math> działania algorytmu w arytmetyce |
| | <math>fl_\nu</math> można zinterpretować jako nieco zaburzony wynik |
| | algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji |
| | <math>y_\nu=(N(f))_\nu\in N(F)</math> o <math>f</math>, przy czym |
| | poziom zaburzeń nie zależy od <math>f</math>. |
|
| |
|
| 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
| | Formalnie znaczy to, że istnieją stałe <math>K_1</math>, <math>K_2</math> oraz |
| działania algorytmu. Jest tak m.in. dla macierzy
| | <math>\nu_0>0</math> takie, że spełniony jest następujący warunek. |
| * symetrycznych, dodatnio określonych: <math>\displaystyle A=A^T</math> oraz <math>\displaystyle x^TAx > 0</math>, <math>\displaystyle \forall x
| | Dla dowolnej <math>\nu\le\nu_0</math> oraz informacji <math>y\in N(F_0)</math> |
| \neq 0</math>, | | można dobrać <math>y_\nu\in N(F)</math> oraz |
| * silnie diagonalnie dominujących: macierz <math>\displaystyle A</math> (lub <math>\displaystyle A^T</math>) spełnia
| | <math>\left(\varphi(y_\nu)\right)_\nu</math> takie, że |
|
| |
|
| <center><math>\displaystyle |a_{ii}| > \sum_{j\neq i} |a_{ij}|, \qquad \forall i. | | <center><math>\|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|, |
| </math></center> | | </math></center> |
|
| |
|
| ==Pamięć hierarchiczna komputerów. Biblioteki BLAS, LAPACK==
| | <center><math>\|\left(\varphi(y_\nu)\right)_\nu - \varphi(y_\nu)\|\,\le\, |
| | K_2\,\nu\,\|\varphi(y_\nu)\|</math>,</center> |
|
| |
|
| W poprzednim wykładzie sprawdziliśmy, jaki jest koszt obliczeniowy algorytmu
| | oraz |
| 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
| | <center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\, |
| większość czasu czekał na dane, a jego efektywność drastycznie spadnie. Z
| | fl_\nu\left(\varphi(N(f))\right)\,=\, |
| niewielką przesadą można powiedzieć, że
| | \left(\varphi(y_\nu)\right)_\nu. |
| | | </math></center> |
| W optymalizacji szybkości działania programu numerycznego,
| |
| obecnie cała walka idzie o to, by procesor przez cały czas ''miał co liczyć''.
| |
|
| |
|
| 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:
| |
|
| |
| ; 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.
| |
|
| |
| 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 korzystający z {cache}'a? Jak go poprawić, by
| |
| korzystał z {cache}'a w sposób wzorowy?===
| |
|
| |
| 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ą {Dongarra} klasyczny już przykład mnożenia dwóch macierzy.
| |
|
| |
| {Kilka implementacji mnożenia dwóch macierzy, a pamięć cache}
| |
|
| |
| 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 !double! 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====
| |
|
| |
| {{algorytm||
| |
| /* 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];
| |
| }} | | }} |
|
| |
|
| Jest to algorytm, który zapewne większości z Czytelników pierwszy przyszedłby
| | [[Image:MNcondition7.png|thumb|550px|center|Numerycznie poprawny algorytm daje w arytmetyce <math>fl_\nu</math> wynik <math>ALG(N(x))</math>, który daje |
| do głowy, gdyż realizuje wprost powszechnie znaną regułę mnożenia macierzy
| | się zinterpretować jako mało zaburzony wynik <math>f(y)</math> zadania na mało zaburzonych |
| "wiersz przez kolumnę". W pamięci {cache} L1 mieści się 64KB danych i
| | danych <math>x</math>.]] |
| 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 !double!.
| |
|
| |
|
| 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>,
| | Zauważmy,że jeśli <math>f\in R^n</math>, |
| {cache miss} następuje praktycznie w każdym kroku. Dzieje się tak dlatego,
| | <math>N(f)=(f_1,\ldots,f_n)</math>, oraz algorytm jest |
| że wymiary naszych macierzy są wielokrotnością 512: odwołując się
| | dokładny, <math>{\bf ALG}\equiv\varphi\equiv S</math>, to numeryczną |
| do kolejnych !B[k*N+j]!, !k! <nowiki>=</nowiki> <math>\displaystyle 0\ldots N</math>, odwołujemy się do co
| | poprawność algorytmu można równoważnie zapisać jako |
| 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====
| | <center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\, |
| | | \left(S(f_\nu)\right)_\nu. |
| Różni się on od poprzedniego jedynie
| |
| kolejnością dwóch wewnętrznych pętli:
| |
| | |
| {{algorytm||
| |
| /* 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];
| |
| }}
| |
| | |
| 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 !B!, 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 !B!. 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>:
| |
| | |
| {{algorytm|| | |
| /* 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];
| |
| | |
| }} | |
| | |
| (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. {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> | | </math></center> |
|
| |
|
| Dla odmiany, C przechowuje w pamięci elementy macierzy wierszami, tzn. kolejno
| | Numeryczna poprawność algorytmu jest wielce pożądaną cechą. |
| <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 !matrix!, gdyż pętle są
| |
| ustawione tak, by przechodzić przez macierz wzdłuż kolejnych kolumn. Dlatego nie
| |
| jest tu konieczne użycie makra !IJ()!, a sprytne wykorzystanie
| |
| pointera !ptr!
| |
| 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 !IJ()! (i obliczania wyrażenia !i+j*N!) 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 !IJ()! 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 {j*N}, 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 {IJ()}) czy też pakiet CLAPACK
| |
| (macierz w wektorze), zob.
| |
| {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|]].
| |
| | |
| ===BLAS, LAPACK i ATLAS===
| |
|
| |
|
| W zadaniach dotyczących macierzy gęstych, warto skorzystać z klasycznego tandemu
| | <blockquote style="background-color: #fefeee; padding:1em; margin-left,margin-right:2em; margin-top,margin-bottom: 1em;"> |
| bibliotek: BLAS ({Basic Linear Algebra Subprograms}) {BLAS-home-page}
| | Algorytm numerycznie poprawny działa w arytmetyce zmiennoprzecinkowej (niemal) tak, jakby wszystkie obliczenia były wykonywane w arytmetyce dokładnej, a tylko dane i wyniki podlegały reprezentacji w skończonej precyzji. |
| oraz LAPACK ({Linear Algebra PACKage}) {LAPACK-home-page}. Dla macierzy
| | </blockquote> |
| rozrzedzonych (zawierających wiele elementów zerowych) warto zastosować bardziej
| |
| wyspecjalizowane biblioteki. Aby wykorzystać do maksimum moc oferowaną przez te
| |
| biblioteki (w połączeniu z naszym komputerem) warto skorzystać z optymalizowanej
| |
| wersji BLASów i LAPACKa, czyli z ATLASa, {ATLAS-home-page}. Istnieje inna
| |
| wersja optymalizowanych BLASów, tzw. Goto BLAS. Niektóre procedury ATLASa są
| |
| istotnie szybsze od standardowych BLASów, dając, przykładowo dla zadania
| |
| mnożenia dwóch macierzy (na komputerze z procesorem Pentium 4 i
| |
| dostatecznie dużych macierzy), '''ponaddziesięciokrotne przyspieszenie''' na
| |
| zmiennych typu {float} i {double} i około pięciokrotne na zmiennych
| |
| typu {complex } i {double complex}.
| |
|
| |
|
| Aby osiągnąć największe przyspieszenie, bibliotekę ATLAS należy skompilować
| | ==Rola uwarunkowania zadania== |
| samemu na własnej (''nieobciążonej'' w trakcie kompilacji!) architekturze. W plikach
| |
| {Makefile} ATLASa brak jednak opcji instalacji bibliotek w standardowych
| |
| miejscach --- trzeba zrobić to samemu.
| |
|
| |
|
| BLAS i LAPACK są często wykorzystywane przez inne biblioteki numeryczne,
| | Niech <math>{\bf ALG}(\cdot)=\varphi(N(\cdot))</math> będzie algorytmem numerycznie |
| na nich opierają się również funkcje macierzowe w Octave i MATLABie.
| | poprawnym dla danych <math>F_0\subset F</math>. Wtedy jego błąd w <math>fl_\nu</math> |
| Optymalizowane biblioteki BLAS i LAPACK dla swoich architektur oferują
| | można oszacować następująco: |
| producenci procesorów Intel (biblioteka MKL) oraz AMD (biblioteka ACML)
| |
|
| |
|
| BLAS {BLAS-home-page} jest kolekcją procedur służących manipulacji podstawowymi obiektami
| | <center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\| \;=\; |
| algebry liniowej: skalarami, wektorami i macierzami. Obsługiwane są obiekty
| | \|S(f)-\left(\varphi(y_\nu)\right)_\nu\| \\ |
| rzeczywiste (w pojedynczej albo podwójnej precyzji) oraz zespolone (także w
| | &\le \|S(f)-\varphi(y)\|\,+\, |
| dwóch precyzjach). W BLAS wyróżnia się 3 poziomy abstrakcji algorytmów:
| | \|\varphi(y)-\varphi(y_\nu)\|\,+\, |
| * BLAS Level 1 -- działania typu wektor-wektor, np. operacja AXPY, czyli
| | \|\varphi(y_\nu)-\left(\varphi(y_\nu)\right)_\nu\| \\ |
| uogólnione
| | &\le \|S(f)-{\bf ALG}(f)\|\,+\, |
| dodawanie wektorów
| | \|\varphi(y)-\varphi(y_\nu)\|\,+\, |
| <center><math>\displaystyle
| | K_2\,\nu\,\|\varphi(y_\nu)\| \\ |
| y \leftarrow \alpha x + y, | | &\le \|S(f)-{\bf ALG}(f)\|\,+\, |
| </math></center> | | (1 + K_2 \nu) \|\varphi(y)-\varphi(y_\nu)\|\,+\, |
| | K_2\,\nu\,\|\varphi(y)\|, |
| | \end{align}</math></center> |
|
| |
|
| albo obliczanie normy wektora. BLAS Level 1 ma za zadanie porządkować kod;
| | przy czym <math>\|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|</math>. Stąd |
| * BLAS Level 2 -- działania typu macierz--wektor, np. mnożenie macierzy
| | w szczególności wynika, że jeśli algorytm jest numerycznie |
| przez wektor
| | poprawny i ciągły ze względu na informację <math>y</math>, to |
| <center><math>\displaystyle
| |
| y \leftarrow \alpha A x + y | |
| </math></center> | |
|
| |
|
| Zapis algorytmów z użyciem BLAS Level 2 umożliwia potencjalnie przyspieszenie
| | <center><math>\lim_{\nu\to 0}\,\|S(f)-fl_\nu({\bf ALG}(f))\|\,=\, |
| programu, m.in. ze względu na to, że zoptymalizowane procedury BLAS Level 2 mogą np.
| | \|S(f)-{\bf ALG}(f)\|. |
| wykorzystywać instrukcje wektorowe nowoczesnych procesorów, zob.
| |
| Rozdział [[##sec:architektura|Uzupelnic sec:architektura|]];
| |
| * BLAS Level 3 -- operacje typu macierz--macierz, np. mnożenie dwóch
| |
| macierzy:
| |
| <center><math>\displaystyle | |
| C \leftarrow \alpha A\cdot B + C
| |
| </math></center> | | </math></center> |
|
| |
|
| W związku z tym, że operacje macierzowe wymagają wykonania <math>\displaystyle O(N^3)</math> działań
| | To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie |
| arytmetycznych przy
| | się zachowywał w <math>fl_\nu</math> prawie tak, jak w arytmetyce idealnej. |
| <math>\displaystyle O(N^2)</math> danych (gdzie <math>\displaystyle N</math> jest wymiarem macierzy), wykorzystanie | |
| zoptymalizowanych procedur BLAS Level 3 może znacząco przyspieszyć wykonywanie
| |
| obliczeń na maszynach z pamięcią hierarchiczną (zob. Rozdział [[##sec:cache|Uzupelnic sec:cache|]]).
| |
|
| |
|
| Podsumowując: przyjmijmy, że dysponujemy dobrze zoptymalizowaną biblioteką BLAS.
| | Z powyższych wzorów wynika, że błąd w <math>fl_\nu</math> algorytmu |
| Wówczas dany algorytm algebry liniowej najlepiej zapisać przy użyciu procedur BLAS Level
| | numerycznie poprawnego zależy w dużym stopniu od: |
| 3, naturalnie, pod warunkiem, że w ogóle daje się to zrobić; typową strategią w
| | * dokładności algorytmu w arytmetyce idealnej, |
| takim wypadku jest tworzenie algorytmów blokowych, operujących na ''blokach''
| | * dokładności <math>\nu</math> arytmetyki <math>fl_\nu</math>, |
| macierzy, zamiast na jej pojedynczych elementach.
| | * wrażliwości algorytmu na małe względne zaburzenia informacji <math>y</math>. |
| | |
| | Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy |
| | trochę więcej uwagi jedynie trzeciemu. |
|
| |
|
| Większość użytkowników BLAS nie będzie jednak miała potrzeby pisania własnych
| | Jeśli <math>\varphi</math> spełnia warunek Lipschitza ze stałą <math>L</math>, |
| algorytmów blokowych, gdyż funkcje rozwiązujące podstawowe zadania numerycznej
| | a dokładniej |
| algebry liniowej: rozwiązywanie układów równań (także nad- i niedookreślonych)
| |
| oraz zadania własnego, znajdują się w doskonałym pakiecie LAPACK
| |
| {LAPACK-home-page}, który
| |
| intensywnie i skutecznie wykorzystuje podprogramy BLAS.
| |
|
| |
|
| Nazwy procedur BLASów i
| | <center><math>\|\varphi(y_\nu)-\varphi(y)\|\,\le\,L\,\|y_\nu-y\|</math>,</center> |
| LAPACKa są cokolwiek enigmatyczne na pierwszy rzut oka, ale w istocie bardzo
| |
| łatwo je odgadywać. Każda nazwa składa się z kilku części, najczęściej jest
| |
| postaci {PRRFF}, gdzie
| |
|
| |
|
| ;
| | to |
| : {P} oznacza precyzję i może przyjmować
| |
| wartości: S,D,C,Z, odpowiadające kolejno pojedynczej i podwójnej precyzji w
| |
| dziedzinie
| |
| rzeczywistej i pojedynczej i podwójnej precyzji w dziedzinie zespolonej;
| |
|
| |
|
| ;
| | <center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\| \\ |
| : {RR} oznacza rodzaj zadania, np. GE oznacza {GEneral}, czyli zadanie ogólne
| | &\le \|S(f)-{\bf ALG}(f)\|\,+\, |
| (praktycznie bez założeń), a SY oznacza {SYmmetric}, czyli zadanie symetryczne; | | (1+K_2\nu)L\|y_\nu-y\|\,+\,K_2\nu\|\varphi(y)\| \\ |
| | &\le \|S(f)-{\bf ALG}(f)\|\,+\, |
| | (1+K_2\nu)LK_1\nu\|y\|\,+\,K_2\nu\|\varphi(y)\|. |
| | \end{align}</math></center> |
|
| |
|
| ;
| | W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny |
| : {FF} wreszcie określa samo zadanie, np. SV oznacza
| | algorytmu proporcjonalnie do <math>\nu</math>. |
| {SolVe} (w domyśle: układ równań), MV --- {Matrix-Vector} (w domyśle: mnożenie),
| |
| EV --- {EigenValues}, czyli wartości własne, itp. Są też warianty
| |
| trzyliterowe, np. TRF ({TRiangular Factorization}) i TRS ({TRiangular
| |
| Solve} --- w domyśle, przy użyciu wcześniej wyznaczonej faktoryzacji)
| |
|
| |
|
| Jeśli jednak ''nie możemy zgadnąć'', jaka jest nazwa procedury BLAS/LAPACKa,
| | Bardziej jednak interesuje nas błąd <strong>względny</strong>. Wybierzmy |
| która byłaby nam potrzebna,
| | "małe" <math>\eta\ge 0</math> i przypuśćmy, że |
| najlepiej przejrzeć (przeszukać) strony internetowe tych pakietów w serwisie
| |
| Netlib.
| |
|
| |
|
| Zestawienie najczęściej wykorzystywanych procedur BLAS i LAPACKa przedstawiamy
| | <center><math>\|\varphi(y_\nu)-\varphi(y)\|\;\le\; |
| poniżej. Każda z tych procedur ma swój wariant "ekspercki", np. dla
| | M\,K_1\,\nu\,\max(\eta,\|\varphi(y)\|)</math>,</center> |
| rozwiązywania układów równań metodą eliminacji Gaussa można skorzystać z
| |
| osobnych procedur generujących rozkład LU oraz z innych, rozwiązujących układy
| |
| trójkątne.
| |
|
| |
|
| {| border=1
| | dla pewnej <math>M</math> niezależnej od <math>y</math>, tzn. błąd względny informacji, |
| |+ <span style="font-variant:small-caps">Uzupelnij tytul</span> | | <math>\|y_\nu-y\|\le K_1\nu\|y\|</math>, przenosi się na błąd względny |
| |-
| | wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia" |
| |
| | <math>M</math>, albo na błąd bezwzględny ze współczynnikiem <math>M\eta</math>. |
| | (Zauważmy, że gdybyśmy wzięli <math>\eta=0</math>, to dla <math>y</math> takiej, że |
| | <math>\varphi(y)=0</math>, musiałoby być <math>\varphi(y_\nu)=0</math> --- co zwykle, choć |
| | nie zawsze, nie jest prawdą.) Wtedy |
|
| |
|
| Zadanie algebry liniowej || Nazwa procedury BLAS/LAPACK
| | <center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\| |
| |- | | & \le \|S(f)-{\bf ALG}(f)\|+ |
| | | | (1 + K_2 \nu) M K_1 \nu \max (\eta, \|\varphi(y)\|)+ |
| | K_2 \nu \|\varphi(y)\| \\ |
| | &= \|S(f)-{\bf ALG}(f)\|\,+\,\nu\, |
| | \Big(\,MK_1(1+K_2\nu)+K_2\Big)\max(\eta,\|\varphi(y)\|). |
| | \end{align}</math></center> |
|
| |
|
| mnożenie wektora przez macierz || DGEMV
| | W szczególności, gdy algorytm jest dokładny i korzysta z pełnej |
| |-
| | informacji o <math>f</math>, tzn. <math>S\equiv{\bf ALG}\equiv\varphi</math>, to |
| | mnożenie macierzy przez macierz || DGEMM
| | błąd |
| |-
| |
| |
| |
| rozwiązywanie układu równań || DGESV
| |
| |-
| |
| | rozkład LU (w miejscu) || DGETRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DGETRF || DGETRS
| |
| |-
| |
| |
| |
| rozwiązywanie układu z macierzą symetryczną || DSYSV
| |
| |-
| |
| | rozkład LDL<math>\displaystyle ^T</math> macierzy symetrycznej (w miejscu) || DSYTRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DSYTRF || DSYTRS
| |
| |-
| |
| |
| |
| rozwiązywanie układu z macierzą pasmową || DGBSV
| |
| |-
| |
| | rozkład LU macierzy pasmowej (w miejscu) || DGBTRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DGBTRF || DGBTRS
| |
| |-
| |
| |
| |
| zagadnienie własne || DGESV
| |
| |-
| |
| |
| |
|
| |
|
| |} | | <center><math>\frac{\|S(f)-fl_\nu({\bf ALG}(f))\|} |
| | | {\max (\eta, \|S(f)\|)} \;\le\; |
| {{przyklad|Mnożenie macierz-wektor w BLAS|| | | \Big( M K_1 (1+K_2\nu) + K_2\Big)\,\nu |
| | | \,\approx\,(M\,K_1\,+\,K_2)\,\nu. |
| Zacznijmy od prostej procedury BLAS Level 2, jaką jest mnożenie macierzy przez
| |
| wektor. Wykorzystamy tzw. wzorcową implementację BLASów (niezoptymalizowaną)
| |
| dostarczaną z dystrybucją np. Red Hat Linux. Jest to biblioteka funkcji
| |
| fortranowskich, dlatego przed dalszym czytaniem warto zapoznać się z
| |
| Rozdziałem [[##sec:FortranC|Uzupelnic sec:FortranC|]].
| |
| | |
| Naszym zadaniem jest wykonanie operacji
| |
| <center><math>\displaystyle
| |
| y \leftarrow \alpha A x + y,
| |
| </math></center> | | </math></center> |
|
| |
|
| gdzie <math>\displaystyle A</math> jest zadaną macierzą <math>\displaystyle N\times M</math>, natomiast <math>\displaystyle y</math> jest wektorem o <math>\displaystyle M</math>
| | Stąd wynika, że jeśli <math>(MK_1+K_2)\nu\ll 1</math>, to błąd względny |
| współrzędnych.
| | algorytmu w <math>fl_\nu</math> jest mały, o ile <math>\|S(f)\|\ge\eta</math>. |
| | Błąd względny jest proporcjonalny do dokładności <math>\nu</math>, |
| | arytmetyki <math>fl_\nu</math>, współczynników proporcjonalności <math>K_i</math> |
| | algorytmu numerycznie poprawnego, oraz do wrażliwości <math>M</math> |
| | zadania <math>S</math> na małe względne zaburzenia danych. |
|
| |
|
| To zadanie realizuje procedura BLASów o nazwie
| | Zwróćmy uwagę na istotny fakt, że interesują nas właściwie |
| {DGEMV}. W rzeczywistości, ta procedura wykonuje ogólniejsze zadanie
| | tylko te zaburzenia danych (informacji), które powstają przy |
| wyznaczania wektora
| | analizie algorytmu numerycznie poprawnego. I tak, jeśli algorytm |
| | jest numerycznie poprawny z pozornymi zaburzeniami danych w normie, |
| | to trzeba zbadać wrażliwość zadania ze względu na zaburzenia |
| | danych w normie. Jeśli zaś mamy pozorne zaburzenia |
| | "po współrzędnych" (co oczywiście implikuje pozorne zaburzenia |
| | w normie) to wystarczy zbadać uwarunkowanie zadania ze względu na |
| | zaburzenia "po współrzędnych", itd. |
|
| |
|
| <center><math>\displaystyle | | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
| y \leftarrow \alpha B x + \beta y,
| | <span style="font-variant:small-caps;">Przykład: Iloczyn skalarny</span> |
| </math></center> | | <div class="solution" style="margin-left,margin-right:3em;"> |
|
| |
|
| przy czym macierz <math>\displaystyle B</math> może być równa albo <math>\displaystyle A</math>, albo <math>\displaystyle A^T</math> (jednak za każdym
| | Załóżmy, że dla danych ciągów rzeczywistych o ustalonej |
| razem argumentem macierzowym, jaki przekazujemy {DGEMV}, jest wyjściowa
| | długości <math>n</math>, <math>a_j</math>, <math>b_j</math>, <math>1\le j\le n</math>, chcemy obliczyć |
| macierz <math>\displaystyle A</math>).
| |
|
| |
|
| Jak wiemy, że jest możliwe bezpośrednie wykorzystanie
| | <center><math>S(a,b)\,=\,\sum_{j=1}^n a_j b_j</math></center> |
| biblioteki fortranowskiej w programie w C, jednak musimy pamiętać, iż macierze z
| |
| jakich ona skorzysta muszą być ułożone ''kolumnami'' w jednolitym bloku
| |
| pamięci.
| |
|
| |
|
| Bazując na opisie procedury {DGEMV} ze
| | Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem |
| strony {opis:dgemv}, w programie w C powinniśmy
| | i korzystający z pełnej informacji o kolejnych współrzednych. |
| napisać prototyp tej funkcji następująco:
| |
|
| |
|
| {{algorytm||
| | Oznaczmy przez <math>\tilde a_j</math> i <math>\tilde b_j</math> reprezentacje liczb |
| int dgemv_( char* TRANS,
| | <math>a_j</math> i <math>b_j</math> w <math>fl_\nu</math>, <math>\tilde a_j=a_j(1+\alpha_j)</math>, |
| int* M,
| | <math>\tilde b_j=b_j(1+\beta_j)</math>, oraz przez <math>\gamma_j</math> i <math>\delta_j</math> |
| int* N,
| | błędy względne powstałe przy kolejnych mnożeniach i dodawaniach. |
| double* ALPHA,
| | Oczywiście <math>|\alpha_j|,|\beta_j|, |\gamma_j|, |\delta_j|\le\nu</math>. |
| double* A,
| | Otrzymujemy |
| int* LDA,
| |
| double* X,
| |
| int* INCX,
| |
| double* BETA,
| |
| double* Y,
| |
| int* INCY );
| |
| }}
| |
| | |
| Dla własnej wygody, a także dla przyszłego wykorzystania, umieścimy ten
| |
| prototyp, razem z innymi przydatnymi drobiazgami (m.in. makro {IJ}
| |
| -- zob. Rozdział [[##sec:macierze-w-komputerze|Uzupelnic sec:macierze-w-komputerze|]] i zmienne całkowite
| |
| !static int BLASONE <nowiki>=</nowiki> 1, BLASMONE <nowiki>=</nowiki> -1;!), w pliku
| |
| nagłówkowym {blaslapack.h}.
| |
| | |
| Dalej już jest łatwo: oto pełny kod programu realizującego operację mnożenia macierzy przez wektor
| |
| przy użyciu procedury BLAS {DGEMV}:
| |
|
| |
|
| {{algorytm|| | | <center><math>\begin{align} fl_\nu\left(\sum_{j=1}^n a_jb_j\right) &= |
| #include <stdio.h>
| | \Big(\,(fl_\nu(\sum_{j=1}^{n-1}a_jb_j)\,+\,\tilde a_n\tilde b_n |
| #include "blaslapack.h"
| | (1+\gamma_n)\,\Big)(1+\delta_n)\,=\,\ldots \\ |
| | &= \bigg(\cdots\Big( |
| | \tilde a_1\tilde b_1(1+\gamma_1)+\tilde a_2\tilde b_2 |
| | (1+\gamma_2)\Big)(1+\delta_2) +\cdots+ |
| | \tilde a_n\tilde b_n(1+\gamma_n)\bigg)(1+\delta_n) \\ |
| | &= \tilde a_1\tilde b_1(1+\gamma_1)(1+\delta_2) |
| | \cdots(1+\delta_n) +\cdots+\tilde a_j |
| | \tilde b_j(1+\gamma_j)(1+\delta_j)\cdots(1+\delta_n) \\ |
| | &= \sum_{j=1}^n a_jb_j(1+e_j), |
| | \end{align}</math></center> |
|
| |
|
| double* mmread(char *filename, int* N, int* M );
| | gdzie w przybliżeniu (tzn. gdy <math>\nu\to 0</math>) mamy <math>|e_1|\leq (n+2)\nu</math> |
| | i <math>|e_j|\leq (n-j+4)\nu</math>, <math>2\le j\le n</math>. Algorytm naturalny jest więc |
| | numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany |
| | w <math>fl_\nu</math> można zinterpretować jako dokładny wynik dla danych |
| | <math>a_{\nu,j}=a_j</math> i <math>b_{\nu,j}=b_j(1+e_j)</math>, przy czym |
| | <math>\|b_\nu-b\|_p\leq (n+2)\nu\|b\|_p</math>. |
|
| |
|
| int main()
| | Zobaczmy teraz, jak błąd we współrzędnych <math>b_j</math> wpływa |
| {
| | na błąd wyniku. Mamy |
| int N, M, i, j;
| |
| double *A, *x, *y;
| |
|
| |
|
| /* wczytujemy macierz z pliku w formacie MatrixMarket */
| | <center><math>\begin{align} \Big|\sum_{j=1}^n a_jb_j-fl_\nu\Big(\sum_{j=1}^n a_jb_j\Big)\Big| |
| /* macierz jest reprezentowana w formacie kolumnowym */
| | &= \Big|\sum_{j=1}^n a_jb_j-\sum_{j=1}^n a_jb_j(1+e_j)\Big| \\ |
| A <nowiki>=</nowiki> mmread("mbeacxc.mtx", &N, &M );
| | &= \Big|\sum_{j=1}^n e_ja_jb_j\Big| |
| | \,\le\, \sum_{j=1}^n |e_j||a_jb_j| \\ |
| | &\leq (n+2)\nu\sum_{j=1}^n |a_jb_j|. |
| | \end{align}</math></center> |
|
| |
|
| x <nowiki>=</nowiki> (double *)malloc(N*sizeof(double));
| | Stąd dla <math>\eta\ge 0</math> |
| y <nowiki>=</nowiki> (double *)malloc(M*sizeof(double));
| |
| for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> N; i++)
| |
| x[IJ(i,1,N)] <nowiki>=</nowiki> (double)i;
| |
|
| |
|
| /* obliczamy y <nowiki>=</nowiki> 5*A*x, korzystając z procedury BLAS L2: DGEMV */
| | <center><math>\frac{|\sum_{j=1}^n a_jb_j-fl_\nu(\sum_{j=1}^n a_jb_j)|} |
| | | {\max(\eta,|\sum_{j=1}^n a_jb_j|)} \,\leq\, |
| { | | K_{\eta}\,(n+2)\,\nu, |
| char TRANS <nowiki>=</nowiki> 'N'; double ALPHA <nowiki>=</nowiki> 5.0, BETA <nowiki>=</nowiki> 0.0;
| |
| | |
| dgemv_(&TRANS, &N, &M, &ALPHA, A, &N, x, &BLASONE,
| |
| &BETA, y, &BLASONE );
| |
| | |
| } | |
| | |
| /* wydruk wyniku */
| |
| for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> M; i++)
| |
| printf("
| |
| return(0);
| |
| } | |
| }}
| |
| | |
| Zwróćmy uwagę na wykorzystanie "stałej" {BLASONE}, równej 1,
| |
| predefiniowanej w pliku {blaslapack.h}. Nasz program kompilujemy
| |
| standardowo, pamiętając o dołączeniu na etapie linkowania używanych przez nas
| |
| bibliotek:
| |
| | |
| gcc -o testblas testblas.c -llapack -lblas -lgfortran -lm
| |
| | |
| --- dokładnie ''w tej właśnie kolejności'' (LAPACK oczywiście w tym momencie
| |
| dołączamy na wyrost: nasz program nie korzysta z żadnej funkcji LAPACKa, wobec
| |
| tego opcja !-llapack! zostanie zignorowana).
| |
| }}
| |
| | |
| Pamiętamy oczywiście, że standardowe BLASy i LAPACK nie są zoptymalizowane w
| |
| stopniu pozwalającym na (prawie) maksymalne wykorzystanie możliwości sprzętu. W
| |
| dzisiejszym świecie, gdzie mamy do czynienia z szerokim wachlarzem procesorów o
| |
| rozmaitych prędkościach, rozmiarach pamięci podręcznej cache i licznymi
| |
| szczegółami architektury komputera, które mają znaczący wpływ na jego
| |
| ostateczne osiągi, wydawać by się mogło, że jedynie dostarczane przez
| |
| producenta BLASy mogłyby osiągnąć zadowalający poziom optymalizacji. Jednakże
| |
| jest inna droga, wykorzystująca kolejne dzieło, jakim jest ATLAS
| |
| ({Automatically Tuned Linear Algebra Subroutines}) {ATLAS-home-page}. ATLAS to faktycznie
| |
| BLAS i (prawie kompletny) LAPACK, a różnica polega na sposobie instalacji tego
| |
| pakietu: otóż w trakcie instalacji automatycznie testuje się sprawność BLAS dla
| |
| różnych wyborów parametrów (np. takich, jak rozmiar bloku macierzy używanego w
| |
| algorytmach blokowych BLAS Level 3) i wybiera ten zestaw parametrów, który
| |
| dawał największą szybkość. Zauważmy więc, że podejście do problemu
| |
| optymalizacji BLAS jest tu dość siłowe (testy mogą trwać nawet kilka godzin!),
| |
| ale dzięki temu, że optymalizacja jest prowadzona drogą eksperymentów, wyniki
| |
| są po prostu znakomite, i to na bardzo wielu różnych architekturach!
| |
| | |
| ====Układy z macierzą ortogonalną====
| |
| | |
| Równie tanio można rozwiązać układ równań
| |
| | |
| <center><math>\displaystyle
| |
| Q\vec x = \vec 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
| |
| \vec x = Q^T \vec 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
| |
| \vec x = Q^H \vec b.
| |
| </math></center> | | </math></center> |
|
| |
|
| ===Metoda eliminacji Gaussa===
| | gdzie |
| | |
| W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym
| |
| rozwiązywania układu równań
| |
| | |
| <center><math>\displaystyle Ax=b</math></center>
| |
|
| |
|
| okazuje się popularna
| | <center><math>K_\eta\,=\,K_\eta(a,b)\,=\,\frac{\sum_{j=1}^n |a_jb_j|} |
| eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm
| | {\max(\eta,|\sum_{j=1}^n a_jb_j|) }. |
| 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> | | </math></center> |
|
| |
|
| a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi: | | Zauważmy, że jeśli iloczyny <math>a_jb_j</math> są wszystkie dodatnie |
| | | albo wszystkie ujemne, to <math>K_\eta=1</math>, tzn. zadanie jest dobrze |
| {{algorytm|Algorytm eliminacji Gaussa|AEG|AEG
| | uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej |
| Znajdź rozkład <math>\displaystyle A=LU</math>;
| | <math>n\nu</math>. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile |
| Rozwiąż <math>\displaystyle Ly = b</math>;
| | liczba <math>n</math> składników nie jest horrendalnie duża. W ogólności |
| Rozwiąż <math>\displaystyle Ux = y</math>;
| | jednak <math>K_\eta</math> może być liczbą dowolnie dużą i wtedy nie możemy |
| }}
| | być pewni uzyskania dobrego wyniku w <math>fl_\nu</math>. |
|
| |
|
| Przypuśćmy, że taki rozkład <math>\displaystyle A=LU</math> istnieje. Wówczas, zapisując macierze w postaci
| | </div></div> |
| blokowej, eskponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy,
| |
| mamy
| |
|
| |
|
| <center><math>\displaystyle | | <!-- |
| \beginpmatrix
| |
| a_{11} & a_{12}^T\\
| |
| a_{21} & A_{22}
| |
| \endpmatrix
| |
| =
| |
| \beginpmatrix
| |
| 1 & 0^T\\
| |
| l_{21} & L_{22}
| |
| \endpmatrix
| |
| \beginpmatrix
| |
| u_{11} & u_{12}^T\\
| |
| 0 & U_{22},
| |
| \endpmatrix
| |
| </math></center>
| |
|
| |
|
| skąd (mnożąc blokowo macierz <math>\displaystyle L</math> przez <math>\displaystyle U</math>) wynika, że
| | <div style="margin-top:1em; padding-top,padding-bottom:1em;"> |
| * <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest
| | <span style="font-variant:small-caps;">Przykład: Pierwiastki trójmianu</span> |
| kopią pierwszego wiersza <math>\displaystyle A</math>,
| | <div class="solution" style="margin-left,margin-right:3em;"> |
| * <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
| | Rozpatrzymy teraz zadanie obliczenia wszystkich pierwiastków |
| rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
| | rzeczywistych równania kwadratowego. |
| klasycznych pętli. Jest to ważne w praktyce numerycznej, gdyż rekurencja
| | Będziemy zakładać, że model obliczeniowy dopuszcza obliczanie |
| kosztuje: zarówno pamięć, jak i czas.
| | pierwiastków kwadratowych z liczb nieujemnych oraz |
| | <math>fl_\nu(\sqrt{x})=rd_\nu(\sqrt{rd_\nu(x)})</math>. |
|
| |
|
| Ponadto zauważmy, że opisany algorytm możemy wykonać w miejscu, nadpisując
| | Okazuje się, że nie umiemy pokazać numerycznej poprawności |
| 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
| | "szkolnego" algorytmu obliczającego pierwiastki równania |
| pamiętać, bo wiemy a priori, że tam są).
| | bezpośrednio ze wzorów omawianych powyżej. Można jednak pokazać |
| | | numeryczną poprawność drobnej jego modyfikacji wykorzystującej |
| {{algorytm|title<nowiki>=</nowiki>Algorytm rozkładu LU||
| | wzory Viete'a. |
| | |
| 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
| |
| }}
| |
| | |
| Ł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ć {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 = \beginpmatrix 0 & 1\\ 1 & 0
| |
| \endpmatrix
| |
| </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 {}{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|title<nowiki>=</nowiki>Algorytm rozkładu LU z wyborem elementu głównego w kolumnie||
| |
| | |
| 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
| |
| }}
| |
| | |
| Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
| |
| już rozkład LU wykonany z wyborem w kolumnie.
| |
| | |
| {{algorytm|title<nowiki>=</nowiki>Algorytm eliminacji Gaussa z wyborem elementu głównego w kolumnie||
| |
| | |
| 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>;
| |
| }}
| |
| | |
| {}{Eliminacja Gaussa z wyborem elementu głównego}
| |
| | |
| 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
| |
| | |
| W optymalizacji szybkości działania programu numerycznego,
| |
| obecnie cała walka idzie o to, by procesor przez cały czas ''miał co liczyć''.
| |
| | |
| 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:
| |
| | |
| ; 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.
| |
| | |
| 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 korzystający z {cache}'a? Jak go poprawić, by
| |
| korzystał z {cache}'a w sposób wzorowy?===
| |
| | |
| 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ą {Dongarra} klasyczny już przykład mnożenia dwóch macierzy.
| |
| | |
| {Kilka implementacji mnożenia dwóch macierzy, a pamięć cache}
| |
| | |
| 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 !double! 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====
| |
| | |
| {{algorytm|||
| |
| /* 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];
| |
| }}
| |
| | |
| 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 !double!.
| |
| | |
| 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 !B[k*N+j]!, !k! <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:
| |
| | |
| {{algorytm|||
| |
| /* 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];
| |
| }}
| |
| | |
| 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 !B!, 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 !B!. 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>:
| |
|
| |
|
| {{algorytm||| | | {{algorytm||| |
| /* bikj(16) */
| | <pre>\EATWSDelta = p*p - q; |
| for (i <nowiki>=</nowiki> 0; i < N; i+<nowiki>=</nowiki>16)
| | if (Delta == 0) |
| for (k <nowiki>=</nowiki> 0; k < N; k+<nowiki>=</nowiki>16)
| | OUT(p); |
| for (j <nowiki>=</nowiki> 0; j < N; j+<nowiki>=</nowiki>16)
| | else |
| for (ii <nowiki>=</nowiki> i; ii < i+15; ii++)
| | if (Delta > 0) |
| for (kk <nowiki>=</nowiki> k; kk < k+15; kk++)
| | { |
| for (jj <nowiki>=</nowiki> j; jj < j+15; jj++)
| | Delta1 = sqrt(d); |
| C[ii*N+jj] +<nowiki>=</nowiki> A[ii*N+kk]*B[kk*N+jj];
| | if (p >= 0) |
| | | { |
| }}
| | x1 = p + Delta1; |
| | | x2 = q/z1; |
| (podobnie działał algorytm bikj(32), tym razem na blokach <math>\displaystyle 32\times 32</math>). | | } |
| | | else |
| 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
| | x2 = p - Delta1; |
| czasie''. Wydzielenie w algorytmie bijk(16) operacji na blokach macierzy ma na celu uniknięcie tego
| | x1 = q/ź2; |
| 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
| | OUT(x1); OUT(x2); |
| algorytmie. Gdyby udało się zachować liczbę {cache misses} na poprzednim poziomie,
| | } |
| można byłoby liczyć na czterokrotne przyspieszenie. {I tak z grubsza jest:
| | </pre>}} |
| 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
| |
| \beginpmatrix
| |
| a_{11} & \cdots & a_{1m}\\
| |
| \vdots & & \vdots\\
| |
| a_{n1} & \cdots & a_{nm}
| |
| \endpmatrix .
| |
| </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 !matrix!, gdyż pętle są
| |
| ustawione tak, by przechodzić przez macierz wzdłuż kolejnych kolumn. Dlatego nie
| |
| jest tu konieczne użycie makra !IJ()!, a sprytne wykorzystanie
| |
| pointera !ptr!
| |
| 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 !IJ()! (i obliczania wyrażenia !i+j*N!) 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 !IJ()! 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 {j*N}, 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 {IJ()}) czy też pakiet CLAPACK
| |
| (macierz w wektorze), zob. | |
| {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|]].
| |
| | |
| ===BLAS, LAPACK i ATLAS===
| |
| | |
| W zadaniach dotyczących macierzy gęstych, warto skorzystać z klasycznego tandemu
| |
| bibliotek: BLAS ({Basic Linear Algebra Subprograms}) {BLAS-home-page}
| |
| oraz LAPACK ({Linear Algebra PACKage}) {LAPACK-home-page}. Dla macierzy
| |
| rozrzedzonych (zawierających wiele elementów zerowych) warto zastosować bardziej
| |
| wyspecjalizowane biblioteki. Aby wykorzystać do maksimum moc oferowaną przez te
| |
| biblioteki (w połączeniu z naszym komputerem) warto skorzystać z optymalizowanej
| |
| wersji BLASów i LAPACKa, czyli z ATLASa, {ATLAS-home-page}. Istnieje inna
| |
| wersja optymalizowanych BLASów, tzw. Goto BLAS. Niektóre procedury ATLASa są
| |
| istotnie szybsze od standardowych BLASów, dając, przykładowo dla zadania
| |
| mnożenia dwóch macierzy (na komputerze z procesorem Pentium 4 i
| |
| dostatecznie dużych macierzy), '''ponaddziesięciokrotne przyspieszenie''' na
| |
| zmiennych typu {float} i {double} i około pięciokrotne na zmiennych
| |
| typu {complex } i {double complex}.
| |
|
| |
|
| Aby osiągnąć największe przyspieszenie, bibliotekę ATLAS należy skompilować
| | Mamy bowiem |
| samemu na własnej (''nieobciążonej'' w trakcie kompilacji!) architekturze. W plikach
| |
| {Makefile} ATLASa brak jednak opcji instalacji bibliotek w standardowych
| |
| miejscach --- trzeba zrobić to samemu.
| |
|
| |
|
| BLAS i LAPACK są często wykorzystywane przez inne biblioteki numeryczne,
| | <center><math>\begin{align} fl_\nu(\Delta(p,q)) &= \Big(p^2(1+\alpha)^2(1+\epsilon_1)-q(1+\beta)\Big) |
| na nich opierają się również funkcje macierzowe w Octave i MATLABie.
| | (1+\epsilon_2) \\ |
| Optymalizowane biblioteki BLAS i LAPACK dla swoich architektur oferują
| | &= \left( p^2-q\frac{(1+\beta)}{(1+\alpha)^2(1+\epsilon_1)}\right) |
| producenci procesorów Intel (biblioteka MKL) oraz AMD (biblioteka ACML)
| | (1+\epsilon_2)(1+\alpha)^2(1+\epsilon_1) \\ |
| | &= \Big(p^2-q(1+\delta)\,\Big)(1+\gamma) \,=\, |
| | \Delta(p,q(1+\delta))(1+\gamma), |
| | \end{align}</math></center> |
|
| |
|
| BLAS {BLAS-home-page} jest kolekcją procedur służących manipulacji podstawowymi obiektami
| | gdzie <math>|\delta|,|\gamma|\leq 4\nu</math>. Wyróżnik obliczony w <math>fl_\nu</math> |
| algebry liniowej: skalarami, wektorami i macierzami. Obsługiwane są obiekty
| | jest więc nieco zaburzonym wyróżnikiem dokładnym dla danych |
| rzeczywiste (w pojedynczej albo podwójnej precyzji) oraz zespolone (także w
| | <math>p</math> i <math>q_\nu=q(1+\delta)</math>. W szczególności |
| dwóch precyzjach). W BLAS wyróżnia się 3 poziomy abstrakcji algorytmów:
| |
| * BLAS Level 1 -- działania typu wektor-wektor, np. operacja AXPY, czyli
| |
| uogólnione
| |
| dodawanie wektorów
| |
| <center><math>\displaystyle | |
| y \leftarrow \alpha x + y,
| |
| </math></center> | |
|
| |
|
| albo obliczanie normy wektora. BLAS Level 1 ma za zadanie porządkować kod;
| | <center><math>\mbox{sgn} (fl_\nu(\Delta(p,q)))= \mbox{sgn} (\Delta(p,q_\nu))</math></center> |
| * BLAS Level 2 -- działania typu macierz--wektor, np. mnożenie macierzy
| |
| przez wektor
| |
| <center><math>\displaystyle | |
| y \leftarrow \alpha A x + y
| |
| </math></center> | |
|
| |
|
| Zapis algorytmów z użyciem BLAS Level 2 umożliwia potencjalnie przyspieszenie
| | Jeśli <math>p\ge 0</math> to |
| programu, m.in. ze względu na to, że zoptymalizowane procedury BLAS Level 2 mogą np.
| |
| wykorzystywać instrukcje wektorowe nowoczesnych procesorów, zob.
| |
| Rozdział [[##sec:architektura|Uzupelnic sec:architektura|]];
| |
| * BLAS Level 3 -- operacje typu macierz--macierz, np. mnożenie dwóch
| |
| macierzy:
| |
| <center><math>\displaystyle
| |
| C \leftarrow \alpha A\cdot B + C
| |
| </math></center> | |
|
| |
|
| W związku z tym, że operacje macierzowe wymagają wykonania <math>\displaystyle O(N^3)</math> działań
| | <center><math>\begin{align} fl_\nu(x1(p,q)) &= \Big(p(1+\alpha)+ |
| arytmetycznych przy
| | \sqrt{fl_\nu(\Delta(p,q))}(1+\epsilon_3)\Big)(1+\epsilon_4) \\ |
| <math>\displaystyle O(N^2)</math> danych (gdzie <math>\displaystyle N</math> jest wymiarem macierzy), wykorzystanie
| | &= \Big(p(1+\alpha)+\sqrt{\Delta(p,q_\nu) |
| zoptymalizowanych procedur BLAS Level 3 może znacząco przyspieszyć wykonywanie
| | (1+\gamma)}(1+\epsilon_3)\Big)(1+\epsilon_4) \\ |
| obliczeń na maszynach z pamięcią hierarchiczną (zob. Rozdział [[##sec:cache|Uzupelnic sec:cache|]]).
| | &= \left(p+\sqrt{\Delta(p,q_\nu)}\frac{\sqrt{1+\gamma}(1+\epsilon_3)} |
| | {1+\alpha}\right)(1+\epsilon_4)(1+\alpha) \\ |
| | &= \left(p+\sqrt{\Delta(p,q_\nu)}\right)(1+e_1), |
| | \end{align}</math></center> |
|
| |
|
| Podsumowując: przyjmijmy, że dysponujemy dobrze zoptymalizowaną biblioteką BLAS.
| | gdzie <math>|e_1|\leq 6\nu</math>. Zauważmy, że ostatnia równość |
| Wówczas dany algorytm algebry liniowej najlepiej zapisać przy użyciu procedur BLAS Level
| | zachodzi dlatego, że dodajemy liczby tego samego znaku. (Inaczej |
| 3, naturalnie, pod warunkiem, że w ogóle daje się to zrobić; typową strategią w
| | <math>|e_1|</math> mogłaby być dowolnie duża i tak byłoby w algorytmie |
| takim wypadku jest tworzenie algorytmów blokowych, operujących na ''blokach''
| | szkolnym.) Dla drugiego pierwiastka mamy |
| macierzy, zamiast na jej pojedynczych elementach.
| |
|
| |
|
| Większość użytkowników BLAS nie będzie jednak miała potrzeby pisania własnych
| | <center><math>fl_\nu(x2(p,q))\,=\,\frac {q(1+\beta)}{fl_\nu(x1(p,q))}(1+\epsilon_5) |
| algorytmów blokowych, gdyż funkcje rozwiązujące podstawowe zadania numerycznej
| | \,=\,\frac{q_\nu}{fl_\nu(x1(p,q))}(1+e_2), |
| algebry liniowej: rozwiązywanie układów równań (także nad- i niedookreślonych)
| |
| oraz zadania własnego, znajdują się w doskonałym pakiecie LAPACK
| |
| {LAPACK-home-page}, który
| |
| intensywnie i skutecznie wykorzystuje podprogramy BLAS.
| |
| | |
| Nazwy procedur BLASów i
| |
| LAPACKa są cokolwiek enigmatyczne na pierwszy rzut oka, ale w istocie bardzo
| |
| łatwo je odgadywać. Każda nazwa składa się z kilku części, najczęściej jest
| |
| postaci {PRRFF}, gdzie
| |
| | |
| ;
| |
| : {P} oznacza precyzję i może przyjmować
| |
| wartości: S,D,C,Z, odpowiadające kolejno pojedynczej i podwójnej precyzji w
| |
| dziedzinie
| |
| rzeczywistej i pojedynczej i podwójnej precyzji w dziedzinie zespolonej;
| |
| | |
| ;
| |
| : {RR} oznacza rodzaj zadania, np. GE oznacza {GEneral}, czyli zadanie ogólne
| |
| (praktycznie bez założeń), a SY oznacza {SYmmetric}, czyli zadanie symetryczne; | |
| | |
| ;
| |
| : {FF} wreszcie określa samo zadanie, np. SV oznacza
| |
| {SolVe} (w domyśle: układ równań), MV --- {Matrix-Vector} (w domyśle: mnożenie),
| |
| EV --- {EigenValues}, czyli wartości własne, itp. Są też warianty
| |
| trzyliterowe, np. TRF ({TRiangular Factorization}) i TRS ({TRiangular
| |
| Solve} --- w domyśle, przy użyciu wcześniej wyznaczonej faktoryzacji)
| |
| | |
| Jeśli jednak ''nie możemy zgadnąć'', jaka jest nazwa procedury BLAS/LAPACKa,
| |
| która byłaby nam potrzebna,
| |
| najlepiej przejrzeć (przeszukać) strony internetowe tych pakietów w serwisie
| |
| Netlib.
| |
| | |
| Zestawienie najczęściej wykorzystywanych procedur BLAS i LAPACKa przedstawiamy
| |
| poniżej. Każda z tych procedur ma swój wariant "ekspercki", np. dla
| |
| rozwiązywania układów równań metodą eliminacji Gaussa można skorzystać z
| |
| osobnych procedur generujących rozkład LU oraz z innych, rozwiązujących układy
| |
| trójkątne.
| |
| | |
| {| border=1
| |
| |+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
| |
| |-
| |
| |
| |
| | |
| Zadanie algebry liniowej || Nazwa procedury BLAS/LAPACK
| |
| |-
| |
| |
| |
| | |
| mnożenie wektora przez macierz || DGEMV
| |
| |-
| |
| | mnożenie macierzy przez macierz || DGEMM
| |
| |-
| |
| |
| |
| rozwiązywanie układu równań || DGESV
| |
| |-
| |
| | rozkład LU (w miejscu) || DGETRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DGETRF || DGETRS
| |
| |-
| |
| |
| |
| rozwiązywanie układu z macierzą symetryczną || DSYSV
| |
| |-
| |
| | rozkład LDL<math>\displaystyle ^T</math> macierzy symetrycznej (w miejscu) || DSYTRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DSYTRF || DSYTRS
| |
| |-
| |
| |
| |
| rozwiązywanie układu z macierzą pasmową || DGBSV
| |
| |-
| |
| | rozkład LU macierzy pasmowej (w miejscu) || DGBTRF
| |
| |-
| |
| | rozwiązywanie układu z rozkładem uzyskanym z DGBTRF || DGBTRS
| |
| |-
| |
| |
| |
| zagadnienie własne || DGESV
| |
| |-
| |
| |
| |
| | |
| |}
| |
| | |
| {{przyklad|Mnożenie macierz-wektor w BLAS||
| |
| | |
| Zacznijmy od prostej procedury BLAS Level 2, jaką jest mnożenie macierzy przez
| |
| wektor. Wykorzystamy tzw. wzorcową implementację BLASów (niezoptymalizowaną)
| |
| dostarczaną z dystrybucją np. Red Hat Linux. Jest to biblioteka funkcji
| |
| fortranowskich, dlatego przed dalszym czytaniem warto zapoznać się z
| |
| Rozdziałem [[##sec:FortranC|Uzupelnic sec:FortranC|]].
| |
| | |
| Naszym zadaniem jest wykonanie operacji
| |
| <center><math>\displaystyle
| |
| y \leftarrow \alpha A x + y,
| |
| </math></center> | | </math></center> |
|
| |
|
| gdzie <math>\displaystyle A</math> jest zadaną macierzą <math>\displaystyle N\times M</math>, natomiast <math>\displaystyle y</math> jest wektorem o <math>\displaystyle M</math> | | gdzie <math>|e_2|\le 8\nu</math>. |
| współrzędnych.
| |
|
| |
|
| To zadanie realizuje procedura BLASów o nazwie
| | Podobny wynik otrzymalibyśmy dla <math>p<0</math>. Algorytm zmodyfikowany |
| {DGEMV}. W rzeczywistości, ta procedura wykonuje ogólniejsze zadanie
| | jest więc numerycznie poprawny, gdyż otrzymane w <math>fl_\nu</math> pierwiastki |
| wyznaczania wektora
| | są nieco zaburzonymi dokładnymi pierwiatkami dla danych |
| | <math>p_\nu=p</math> i <math>q_\nu=q(1+\delta)</math>. |
|
| |
|
| <center><math>\displaystyle | | Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie |
| y \leftarrow \alpha B x + \beta y,
| | zadania ze względu na zaburzenie danej <math>q</math>, ponieważ pokazaliśmy, |
| </math></center> | | że zaburzenia <math>p</math> można przenieść na zaburzenia <math>q</math> i wyniku. |
| | Niestety, choć algorytm jest numerycznie poprawny, zaburzenia |
| | <math>q</math> mogą sprawić, że nawet znak wyróżnika <math>\Delta</math> może być |
| | obliczony nieprawidłowo. Na przykład dla <math>p=1</math> i <math>q=1\pm 10^{t+1}</math> |
| | mamy <math>\Delta(p,q)=\mp 10^{t+1}</math>, ale |
| | <math>\Delta(rd_\nu(p),rd_\nu(q))=\Delta(1,1)=0</math>. Ogólnie |
|
| |
|
| przy czym macierz <math>\displaystyle B</math> może być równa albo <math>\displaystyle A</math>, albo <math>\displaystyle A^T</math> (jednak za każdym
| | <center><math>|fl_\nu(\Delta(p,q))-\Delta(p,q)|\,\leq\,4\nu(p^2+2|q|)</math>,</center> |
| razem argumentem macierzowym, jaki przekazujemy {DGEMV}, jest wyjściowa
| |
| macierz <math>\displaystyle A</math>).
| |
|
| |
|
| Jak wiemy, że jest możliwe bezpośrednie wykorzystanie
| | a więc tylko dla <math>|\Delta(p,q)|=|p^2-q|>4\nu (p^2+2|q|)</math> |
| biblioteki fortranowskiej w programie w C, jednak musimy pamiętać, iż macierze z
| | możemy być pewni obliczenia właściwego znaku <math>\Delta</math>. Przy |
| jakich ona skorzysta muszą być ułożone ''kolumnami'' w jednolitym bloku
| | tym warunku oraz <math>\Delta>0</math> błąd danych przenosi się w |
| pamięci.
| | normie euklidesowej na błąd wyniku następująco: |
|
| |
|
| Bazując na opisie procedury {DGEMV} ze
| | <center><math>\begin{align} \lefteqn{ \Big( (x1(p,q) - x1(p,q_\nu))^2 |
| strony {opis:dgemv}, w programie w C powinniśmy
| | +(x2(p,q) - x2(p,q_\nu))^2 \Big)^{1/2} } \\ |
| napisać prototyp tej funkcji następująco:
| | &= \frac{\sqrt 2 |\delta q|} {\sqrt{p^2-q}+\sqrt{p^2-q_\nu}} |
| | \,\leq\, 2\sqrt 2 \nu\frac{|q|}{\sqrt{p^2-q}} \\ |
| | &= 2\sqrt 2 \nu\frac{|q|/p^2}{\sqrt{1-q/p^2} |
| | \max(\eta/|p|,\sqrt{2(1+(1-q/p^2))}) } \\ |
| | & & \qquad\qquad\qquad\cdot\max(\eta,(x1(p,q)^2+x2(p,q)^2)^{1/2}). |
| | \end{align}</math></center> |
|
| |
|
| {{algorytm|||
| | Stąd widać, że zadanie jest dobrze uwarunkowane dla <math>q/p^2\ll 1</math> |
| int dgemv_( char* TRANS,
| | i może być źle uwarunkowane dla <math>q/p^2\approx 1</math>. W ostatnim |
| int* M,
| | przypadku nie możemy być pewni otrzymania dobrego wyniku w <math>fl_\nu</math>. |
| int* N,
| | </div></div> |
| double* ALPHA,
| |
| double* A,
| |
| int* LDA,
| |
| double* X,
| |
| int* INCX,
| |
| double* BETA,
| |
| double* Y,
| |
| int* INCY );
| |
| }}
| |
| | |
| Dla własnej wygody, a także dla przyszłego wykorzystania, umieścimy ten
| |
| prototyp, razem z innymi przydatnymi drobiazgami (m.in. makro {IJ}
| |
| -- zob. Rozdział [[##sec:macierze-w-komputerze|Uzupelnic sec:macierze-w-komputerze|]] i zmienne całkowite
| |
| !static int BLASONE <nowiki>=</nowiki> 1, BLASMONE <nowiki>=</nowiki> -1;!), w pliku
| |
| nagłówkowym {blaslapack.h}.
| |
| | |
| Dalej już jest łatwo: oto pełny kod programu realizującego operację mnożenia macierzy przez wektor
| |
| przy użyciu procedury BLAS {DGEMV}:
| |
| | |
| {{algorytm|||
| |
| <nowiki>
| |
| #include <stdio.h>
| |
| #include "blaslapack.h"
| |
| | |
| double* mmread(char *filename, int* N, int* M );
| |
| | |
| int main()
| |
| {
| |
| int N, M, i, j;
| |
| double *A, *x, *y;
| |
| | |
| /* wczytujemy macierz z pliku w formacie MatrixMarket */
| |
| /* macierz jest reprezentowana w formacie kolumnowym */
| |
| A <nowiki>=</nowiki> mmread("mbeacxc.mtx", &N, &M );
| |
| | |
| x <nowiki>=</nowiki> (double *)malloc(N*sizeof(double));
| |
| y <nowiki>=</nowiki> (double *)malloc(M*sizeof(double));
| |
| for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> N; i++)
| |
| x[IJ(i,1,N)] <nowiki>=</nowiki> (double)i;
| |
| | |
| /* obliczamy y <nowiki>=</nowiki> 5*A*x, korzystając z procedury BLAS L2: DGEMV */
| |
| | |
| {
| |
| char TRANS <nowiki>=</nowiki> 'N'; double ALPHA <nowiki>=</nowiki> 5.0, BETA <nowiki>=</nowiki> 0.0;
| |
| | |
| dgemv_(&TRANS, &N, &M, &ALPHA, A, &N, x, &BLASONE,
| |
| &BETA, y, &BLASONE );
| |
| | |
| }
| |
|
| |
|
| /* wydruk wyniku */
| | --> |
| for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> M; i++)
| |
| printf("
| |
| return(0);
| |
| }
| |
| }}
| |
|
| |
|
| Zwróćmy uwagę na wykorzystanie "stałej" {BLASONE}, równej 1,
| | ==Literatura== |
| predefiniowanej w pliku {blaslapack.h}. Nasz program kompilujemy
| |
| standardowo, pamiętając o dołączeniu na etapie linkowania używanych przez nas
| |
| bibliotek:
| |
|
| |
|
| gcc -o testblas testblas.c -llapack -lblas -lgfortran -lm
| | W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 2.3</b> w |
| | * D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X. |
|
| |
|
| --- dokładnie ''w tej właśnie kolejności'' (LAPACK oczywiście w tym momencie
| | Warto także przejrzeć rozdział 2 w |
| dołączamy na wyrost: nasz program nie korzysta z żadnej funkcji LAPACKa, wobec
| | * <span style="font-variant:small-caps">P. Deulfhard, A. Hohmann</span>, <cite>Numerical Analysis in Modern Scientific Computing</cite>, Springer, 2003, |
| tego opcja !-llapack! zostanie zignorowana).
| |
| }}
| |
|
| |
|
| Pamiętamy oczywiście, że standardowe BLASy i LAPACK nie są zoptymalizowane w
| | omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów. |
| stopniu pozwalającym na (prawie) maksymalne wykorzystanie możliwości sprzętu. W
| | Nieocenioną monografią na ten temat jest |
| dzisiejszym świecie, gdzie mamy do czynienia z szerokim wachlarzem procesorów o
| | * <span style="font-variant:small-caps">N. Higham</span>, <cite>Accuracy and Stability of Numerical Algorithms</cite>, SIAM, 2002. |
| rozmaitych prędkościach, rozmiarach pamięci podręcznej cache i licznymi
| |
| szczegółami architektury komputera, które mają znaczący wpływ na jego
| |
| ostateczne osiągi, wydawać by się mogło, że jedynie dostarczane przez
| |
| producenta BLASy mogłyby osiągnąć zadowalający poziom optymalizacji. Jednakże
| |
| jest inna droga, wykorzystująca kolejne dzieło, jakim jest ATLAS | |
| ({Automatically Tuned Linear Algebra Subroutines}) {ATLAS-home-page}. ATLAS to faktycznie
| |
| BLAS i (prawie kompletny) LAPACK, a różnica polega na sposobie instalacji tego
| |
| pakietu: otóż w trakcie instalacji automatycznie testuje się sprawność BLAS dla
| |
| różnych wyborów parametrów (np. takich, jak rozmiar bloku macierzy używanego w
| |
| algorytmach blokowych BLAS Level 3) i wybiera ten zestaw parametrów, który
| |
| dawał największą szybkość. Zauważmy więc, że podejście do problemu
| |
| optymalizacji BLAS jest tu dość siłowe (testy mogą trwać nawet kilka godzin!),
| |
| ale dzięki temu, że optymalizacja jest prowadzona drogą eksperymentów, wyniki
| |
| są po prostu znakomite, i to na bardzo wielu różnych architekturach!
| |