MN04: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
m Zastępowanie tekstu – „,↵</math>” na „</math>,”
 
(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/&nbsp;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ł&nbsp;[[##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ł&nbsp;[[##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&nbsp;MHz <nowiki>=</nowiki>  100&nbsp;MHz/10 (bo
magistrala ({bus}) jest taktowana 100&nbsp;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&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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&nbsp;MHz <nowiki>=</nowiki>  100&nbsp;MHz/10 (bo
magistrala ({bus}) jest taktowana 100&nbsp;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&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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ł&nbsp;[[##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&nbsp;[[##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ł&nbsp;[[##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!

Aktualna wersja na dzień 21:46, 11 wrz 2023


Własności zadania obliczeniowego i algorytmu numerycznego

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Uwarunkowanie zadania obliczeniowego

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 1/10)
  • błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie f(x)=a, ale a 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 wpływu zaburzenia danych na wynik jest kapitalna w świecie obliczeń numerycznych w ogólności, a w szególności --- inżynierskich.

Wprowadza się pojęcie uwarunkowania zadania, to znaczy jego podatności na zaburzenia danych. Dla przejrzystości przypuśćmy, że nasze zadanie obliczeniowe polega na wyznaczeniu f(x) dla danego x.


Jak bardzo będzie odległe f(x~), gdy x~x? Rozważa się dwa przypadki:

  • uwarunkowanie względne: jak względne zaburzenie danych wpływa na błąd względny wyniku:
    ||f(x)f(x~)||||f(x)||condrel(f,x)||xx~||||x||

Najmniejszy mnożnik condrel(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia f(x) dla danego x.

  • uwarunkowanie bezwzględne: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku:
    ||f(x)f(x~)||condabs(f,x)||xx~||

Najmniejszy mnożnik condabs(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia f(x) dla danego x.

Powiemy, że zadanie f(x) jest

  • dobrze uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle postawione w punkcie x, gdy cond(f,x)=+.

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!

Przykład: Uwarunkowanie zadania obliczenia sumy

Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia s(x,y)=x+y ma

condabs(s,(a,b))=1,condrel(s,(a,b))=|a|+|b||a+b|

Tak więc, gdy ab, to condrel(s,(a,b))+ i zadanie jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego, najczęściej rzeczywiście tak będzie...

Przykład

Dla różniczkowalnej funkcji skalarnej f:RR mamy

|f(x)f(x~)||f(x)||xx~|

i w konsekwencji dla zadania obliczenia f(x) dla danego x mamy, przy założeniu małych zaburzeń,

condabs(f,x)=|f(x)|,condrel(f,x)=|f(x)||x||f(x)|

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 złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, ale wręcz pomaga szybciej rozwiązać zadanie główne!

Rozkład algorytmu względem informacji

Algorytm to dokładnie określona i dozwolona w danym modelu obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego zadania (w sposób dokładny lub przybliżony).

Z każdym algorytmem związany jest operator

𝐀𝐋𝐆:FG,

taki że 𝐀𝐋𝐆(f) jest wynikiem działania algorytmu w arytmetyce idealnej dla danej f.

Zauważmy, że wynik 𝐀𝐋𝐆(f) działania algorytmu nie zależy bezpośrednio od f, ale raczej od informacji o f (uzyskanej dzięki poleceniu 𝒩). Informacja ta może być pełna albo tylko częściowa. Informacja jest pełna gdy, np. f=(f1,,fn)Rn i wczytamy wszystkie współrzędne fi. Informacja może być częściowa, gdy f jest funkcją. Wtedy wiele danych może posiadać tę samą informację, co łatwo zaobserwować na przykładzie zadania całkowania.

Niech N:Fn=0Rn będzie operatorem informacji, tzn.

N(f)=(y1,y2,,yn)

jest informacją o f zebraną przy idealnej realizacji algorytmu. Zauważmy, że informacja jest pełna, gdy N jest przekształceniem różnowartościowym, tzn. jeśli f1f2 implikuje N(f1)N(f2). W przeciwnym przypadku mamy do czynienia z informacją częściową.

Każdy algorytm 𝐀𝐋𝐆 może być przedstawiony jako złożenie operatora informacji i pewnego operatora φ:N(F)G zdefiniowanego równością

φ(N(f))=𝐀𝐋𝐆(f).

Zauważmy, że w przypadku informacji częściowej zwykle nie istnieje algorytm dający dokładne rozwiązanie zadania dla każdej danej fF, ponieważ dla danych o tej samej informacji mogą istnieć różne rozwiązania.

Problem wyboru algorytmu

Wybór algorytmu jest najistotniejszą częścią całego procesu numerycznego rozwiązywania zadania. Kierujemy się przy tym przede wszystkim następującymi kryteriami:

  • dokładnością algorytmu,
  • złożonością algorytmu,
  • własnościami numerycznymi algorytmu.

Przez dokładność algorytmu rozumiemy różnicę między rozwiązaniem dokładnym S(f) a rozwiązaniem 𝐀𝐋𝐆(f) dawanym przez algorytm w arytmetyce idealnej. Jeśli 𝐀𝐋𝐆(f)=S(f), fF, to algorytm nazywamy dokładnym.

Mówiąc o złożoności, mamy na myśli złożoność pamięciową (zwykle jest to liczba stałych i zmiennych używanych przez algorytm), jak również złożoność obliczeniową. Na złożoność obliczeniową algorytmu dla danej f składa się koszt uzyskania infomacji y=N(f) (zwykle jest on proporcjonalny do liczby wywołań polecenia 𝒩), oraz koszt kombinatoryczny przetworzenia tej informacji, aż do uzyskania wyniku φ(y). Koszt kombinatoryczny zwykle mierzymy liczbą operacji arytmetycznych wykonywanych przez algorytm.

Przez własności numeryczne algorytmu rozumiemy jego własności przy realizacji w arytmetyce flν. Temu ważnemu tematowi poświęcimy teraz osobny paragraf.

Numeryczna poprawność algorytmu

Pożądane jest, aby algorytm dawał "dobry" wynik zarówno w arytmetyce idealnej, jak i w arytmetyce flν. Niestety, jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm jest dokładny, to w wyniku jego realizacji w flν możemy otrzymać wynik flν(𝐀𝐋𝐆(f)) daleko odbiegający od S(f). W szczególności, prawie zawsze mamy

S(f)flν(𝐀𝐋𝐆(f)).

Zauważmy również, że o ile z reguły znamy dokładne zachowanie się algorytmu w arytmetyce idealnej dla danej informacji, to nie można tego samego powiedzieć o jego zachowaniu się w arytmetyce flν. W związku z tym powstaje pytanie, jak kontrolować błąd algorytmów wynikający z błędów zaokrągleń i jakie algorytmy uznamy za te o najwyższej jakości numerycznej.

Istnienie błędów reprezentacji liczb rzeczywistych powoduje, że informacja y=N(f) o danej f nie jest w ogólności reprezentowana dokładnie. Znaczy to, że zamiast na informacji dokładnej, dowolny algorytm będzie operować na informacji nieco zaburzonej yν, tzn. zaburzonej na poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm będzie w ogólności zaburzony na poziomie błędu reprezentacji. W najlepszym więc wypadku wynikiem działania algorytmu w flν będzie (φ(yν))ν zamiast φ(y). Algorytmy dające tego rodzaju wyniki uznamy za posiadające najlepsze własności numeryczne w arytmetyce flν i nazwiemy numerycznie poprawnymi.

Powiemy, że ciąg rzeczywisty aν=(aν,1,,aν,n) (a właściwie rodzina ciągów {aν}ν) jest nieco zaburzonym ciągiem a=(a1,,an), jeśli istnieje stała K taka, że dla wszystkich dostatecznie małych ν zachodzi

|aν,jaj|Kν|aj|,1jn,

albo ogólniej

aνaKνa,

gdzie jest pewną normą w Rn. W pierwszym przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim o zaburzeniu w normie .

Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas

aνa=max1jn|aν,jaj|Kνmax1jn|aj|=Kνa,

i korzystając z faktu, że w przestrzeni skończenie wymiarowej wszystkie normy są równoważne, otrzymujemy dla pewnych stałych K1 i K2

aνaK1aνaK1KνaK2K1Kνa,

czyli nierówność dla zaburzenia w normie, ze stałą K=K2K1K.

Definicja Algorytm numerycznie poprawny

Algorytm 𝐀𝐋𝐆 rozwiązywania zadania nazywamy numerycznie poprawnym w zbiorze danych F0F, jeśli dla każdej danej fF0 wynik flν(𝐀𝐋𝐆(f)) działania algorytmu w arytmetyce flν można zinterpretować jako nieco zaburzony wynik algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji yν=(N(f))νN(F) o f, przy czym poziom zaburzeń nie zależy od f.

Formalnie znaczy to, że istnieją stałe K1, K2 oraz ν0>0 takie, że spełniony jest następujący warunek. Dla dowolnej νν0 oraz informacji yN(F0) można dobrać yνN(F) oraz (φ(yν))ν takie, że

yνyK1νy,
(φ(yν))νφ(yν)K2νφ(yν),

oraz

flν(𝐀𝐋𝐆(f))=flν(φ(N(f)))=(φ(yν))ν.
Numerycznie poprawny algorytm daje w arytmetyce flν wynik ALG(N(x)), który daje się zinterpretować jako mało zaburzony wynik f(y) zadania na mało zaburzonych danych x.

Zauważmy,że jeśli fRn, N(f)=(f1,,fn), oraz algorytm jest dokładny, 𝐀𝐋𝐆φS, to numeryczną poprawność algorytmu można równoważnie zapisać jako

flν(𝐀𝐋𝐆(f))=(S(fν))ν.

Numeryczna poprawność algorytmu jest wielce pożądaną cechą.

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.

Rola uwarunkowania zadania

Niech 𝐀𝐋𝐆()=φ(N()) będzie algorytmem numerycznie poprawnym dla danych F0F. Wtedy jego błąd w flν można oszacować następująco:

S(f)flν(𝐀𝐋𝐆(f))=S(f)(φ(yν))νS(f)φ(y)+φ(y)φ(yν)+φ(yν)(φ(yν))νS(f)𝐀𝐋𝐆(f)+φ(y)φ(yν)+K2νφ(yν)S(f)𝐀𝐋𝐆(f)+(1+K2ν)φ(y)φ(yν)+K2νφ(y),

przy czym yνyK1νy. Stąd w szczególności wynika, że jeśli algorytm jest numerycznie poprawny i ciągły ze względu na informację y, to

limν0S(f)flν(𝐀𝐋𝐆(f))=S(f)𝐀𝐋𝐆(f).

To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie się zachowywał w flν prawie tak, jak w arytmetyce idealnej.

Z powyższych wzorów wynika, że błąd w flν algorytmu numerycznie poprawnego zależy w dużym stopniu od:

  • dokładności algorytmu w arytmetyce idealnej,
  • dokładności ν arytmetyki flν,
  • wrażliwości algorytmu na małe względne zaburzenia informacji y.

Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy trochę więcej uwagi jedynie trzeciemu.

Jeśli φ spełnia warunek Lipschitza ze stałą L, a dokładniej

φ(yν)φ(y)Lyνy,

to

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)Lyνy+K2νφ(y)S(f)𝐀𝐋𝐆(f)+(1+K2ν)LK1νy+K2νφ(y).

W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny algorytmu proporcjonalnie do ν.

Bardziej jednak interesuje nas błąd względny. Wybierzmy "małe" η0 i przypuśćmy, że

φ(yν)φ(y)MK1νmax(η,φ(y)),

dla pewnej M niezależnej od y, tzn. błąd względny informacji, yνyK1νy, przenosi się na błąd względny wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia" M, albo na błąd bezwzględny ze współczynnikiem Mη. (Zauważmy, że gdybyśmy wzięli η=0, to dla y takiej, że φ(y)=0, musiałoby być φ(yν)=0 --- co zwykle, choć nie zawsze, nie jest prawdą.) Wtedy

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)MK1νmax(η,φ(y))+K2νφ(y)=S(f)𝐀𝐋𝐆(f)+ν(MK1(1+K2ν)+K2)max(η,φ(y)).

W szczególności, gdy algorytm jest dokładny i korzysta z pełnej informacji o f, tzn. S𝐀𝐋𝐆φ, to błąd

S(f)flν(𝐀𝐋𝐆(f))max(η,S(f))(MK1(1+K2ν)+K2)ν(MK1+K2)ν.

Stąd wynika, że jeśli (MK1+K2)ν1, to błąd względny algorytmu w flν jest mały, o ile S(f)η. Błąd względny jest proporcjonalny do dokładności ν, arytmetyki flν, współczynników proporcjonalności Ki algorytmu numerycznie poprawnego, oraz do wrażliwości M zadania S na małe względne zaburzenia danych.

Zwróćmy uwagę na istotny fakt, że interesują nas właściwie tylko te zaburzenia danych (informacji), które powstają przy 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.

Przykład: Iloczyn skalarny

Załóżmy, że dla danych ciągów rzeczywistych o ustalonej długości n, aj, bj, 1jn, chcemy obliczyć

S(a,b)=j=1najbj

Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem i korzystający z pełnej informacji o kolejnych współrzednych.

Oznaczmy przez a~j i b~j reprezentacje liczb aj i bj w flν, a~j=aj(1+αj), b~j=bj(1+βj), oraz przez γj i δj błędy względne powstałe przy kolejnych mnożeniach i dodawaniach. Oczywiście |αj|,|βj|,|γj|,|δj|ν. Otrzymujemy

flν(j=1najbj)=((flν(j=1n1ajbj)+a~nb~n(1+γn))(1+δn)==((a~1b~1(1+γ1)+a~2b~2(1+γ2))(1+δ2)++a~nb~n(1+γn))(1+δn)=a~1b~1(1+γ1)(1+δ2)(1+δn)++a~jb~j(1+γj)(1+δj)(1+δn)=j=1najbj(1+ej),

gdzie w przybliżeniu (tzn. gdy ν0) mamy |e1|(n+2)ν i |ej|(nj+4)ν, 2jn. Algorytm naturalny jest więc numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany w flν można zinterpretować jako dokładny wynik dla danych aν,j=aj i bν,j=bj(1+ej), przy czym bνbp(n+2)νbp.

Zobaczmy teraz, jak błąd we współrzędnych bj wpływa na błąd wyniku. Mamy

|j=1najbjflν(j=1najbj)|=|j=1najbjj=1najbj(1+ej)|=|j=1nejajbj|j=1n|ej||ajbj|(n+2)νj=1n|ajbj|.

Stąd dla η0

|j=1najbjflν(j=1najbj)|max(η,|j=1najbj|)Kη(n+2)ν,

gdzie

Kη=Kη(a,b)=j=1n|ajbj|max(η,|j=1najbj|).

Zauważmy, że jeśli iloczyny ajbj są wszystkie dodatnie albo wszystkie ujemne, to Kη=1, tzn. zadanie jest dobrze uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej nν. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile liczba n składników nie jest horrendalnie duża. W ogólności jednak Kη może być liczbą dowolnie dużą i wtedy nie możemy być pewni uzyskania dobrego wyniku w flν.


Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 2.3 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Warto także przejrzeć rozdział 2 w

  • P. Deulfhard, A. Hohmann, Numerical Analysis in Modern Scientific Computing, Springer, 2003,

omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów. Nieocenioną monografią na ten temat jest

  • N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002.