MN04: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
← poprzednia edycjanastępna edycja →
WizualnieWikikod
Przykry (dyskusja | edycje)
647 edycji
mNie podano opisu zmian
Przykry (dyskusja | edycje)
647 edycji
mNie podano opisu zmian
Linia 2: Linia 2:
''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki''
''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki''


==Rozwiązywanie układów równań liniowych==
{{algorytm|||
 
UULLLULA
Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń
<pre>
zmiennoprzecinkowych, przy czym wiele znanych, ''matematycznie równoważnych'' metod
rozwiązywania takich zadań, ma ''diametralnie różne'' własności numeryczne.
Bardzo ważną klasą takich zadań jest rozwiązywanie układów równań liniowych
 
<center><math>\displaystyle
Ax = b,
</math></center>
 
gdzie <math>\displaystyle A</math> jest nieosobliwą macierzą <math>\displaystyle N\times N</math>, a dany wektor prawej strony <math>\displaystyle b\in R^N</math>.
 
W
praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
 
Rozwiązywanie układów równań liniowych jest sercem wielu innych algorytmów
numerycznych, dlatego nie dziwi, że szacuje się, że około 75 procent czasu
obliczeniowego superkomputerów jest wykorzystywanych właśnie na rozwiązywanie
takich zadań.
 
Okazuje się, że kilka znanych w matematyce sposobów rozwiązywania układów równań
liniowych, takie jak:
* metoda wyznacznikowa (wzory Cramera)
* obliczenie macierzy <math>\displaystyle A^{-1}</math> i następnie <math>\displaystyle x = A^{-1}b</math>
'''nie nadają się''' do numerycznego rozwiązywania takich zadań.
 
O tym, jak ''skutecznie'' rozwiązywać takie zadania, jakie jest ich
uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z
dobrych metod jest metoda eliminacji Gaussa.
 
===Proste układy równań===
 
Niektóre układy równań można bardzo łatwo rozwiązać. Zgodnie z zasadą
 
<blockquote> 
Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
</blockquote>
 
w dalszej kolejności pokażemy, jak dowolny układ równań sprowadzić do sekwencji
dwóch (czasem trzech) łatwych do rozwiązania układów równań. Ale... jakie układy
równań są "łatwe"?
 
====Układy z macierzą trójkątną====
 
Rozważmy układ z macierzą
trójkątną <math>\displaystyle A</math>. Będą nas szczególnie interesować macierze
''trójkątne górne'', dla których <math>\displaystyle a_{i,j}=0</math> gdy <math>\displaystyle i>j</math>, oraz
macierze ''trójkątne dolne'' z jedynkami na przekątnej, tzn.
<math>\displaystyle a_{i,j}=0</math>, <math>\displaystyle i<j</math>, oraz <math>\displaystyle a_{i,i}=1</math>. Macierze pierwszego rodzaju
będziemy oznaczać przez <math>\displaystyle U</math>, a drugiego rodzaju przez <math>\displaystyle L</math>.
 
<center><math>\displaystyle L = \begin{pmatrix}
1 &  &  &        &  &  \\
* & 1 &  &        &  &  \\
* & * & 1 &        &  &  \\
* & * & * & 1 &  &        \\
\vdots  & \vdots  & \vdots  & \ddots  & \ddots &  \\
*  &  *  & * &  \cdots  &  *    & 1
\end{pmatrix} ,
\qquad
U = \begin{pmatrix}
* & * & * & *      & \cdots & * \\
  & * & * & *      & \cdots & * \\
  &  & * & *      & \cdots & * \\
  &  &        & * & \ddots &  \vdots \\
  &  &  &        & \ddots & * \\
  &  &  &        &        & * \end{pmatrix}
</math></center>
 
Układ z nieosobliwą macierzą trójkątną górną
 
<center><math>\displaystyle
  U\, x\;=\; c,
</math></center>
 
<math>\displaystyle U=(u_{i,j})</math>, <math>\displaystyle  c=(c_j)</math>, można rozwiązać stosując algorytm:
<pre>  
<math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
<math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
for (i <nowiki>=</nowiki> N-1; i ><nowiki>=</nowiki> 1; i--)
for (i <nowiki>=</nowiki> N-1; i ><nowiki>=</nowiki> 1; i--)
<math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
<math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
u_{i,j}x_j^*\right)/u_{i,i}</math>;
u_{i,j}x_j^*\right)/u_{i,i}</math>;
</pre>  
</pre>}}
 
(Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy
implikuje, że <math>\displaystyle u_{i,i}\ne 0</math>, <math>\displaystyle \forall i</math>.) Podobnie, układ
<math>\displaystyle L x= c</math> rozwiązujemy algorytmem:
 
<pre>
<math>\displaystyle x_1 = c_1</math>;
for (i<nowiki>=</nowiki>2; i <<nowiki>=</nowiki> N; i++)
<math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;
</pre>
 
Oba algorytmy wymagają rzędu <math>\displaystyle N^2/2</math> mnożeń lub dzieleń i
<math>\displaystyle N^2/2</math> dodawań lub odejmowań, a więc łącznie <math>\displaystyle O(N^2)</math> 
działań arytmetycznych.
 
====Układy z macierzą ortogonalną====
 
Równie tanio można rozwiązać układ równań
 
<center><math>\displaystyle
Q x =  b,
</math></center>
 
gdy <math>\displaystyle Q</math> jest macierzą ortogonalną, to znaczy <math>\displaystyle Q^TQ = I</math>. Rzeczywiście, z
ortogonalności mamy natychmiast, że
 
<center><math>\displaystyle
x = Q^T  b
</math></center>
 
i w konsekwencji <math>\displaystyle x</math> można wyznaczyć kosztem takim jak koszt mnożenia macierzy
przez wektor, czyli <math>\displaystyle O(N^2)</math> operacji.
 
Podobnie, gdy <math>\displaystyle Q\in C^{N\times N}</math> jest unitarna, to znaczy <math>\displaystyle Q^HQ = I</math>,
rozwiązaniem układu równań jest
 
<center><math>\displaystyle
x = Q^H  b.
</math></center>
 
===Metoda eliminacji Gaussa===
 
W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym
rozwiązywania układu równań
 
<center><math>\displaystyle Ax=b</math></center>
 
okazuje się popularna
eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm
ten wyrazimy w języku tzw. ''rozkładu LU'' macierzy, to znaczy,
sprowadzającego zadanie do znalezienia
macierzy trójkątnej dolnej <math>\displaystyle L</math> (z jedynkami na diagonali) oraz trójkątnej górnej
<math>\displaystyle U</math> takich, że
<center><math>\displaystyle
A = LU,
</math></center>
 
a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:
 
<pre> [Algorytm eliminacji Gaussa]
Znajdź rozkład <math>\displaystyle A=LU</math>;
Rozwiąż <math>\displaystyle Ly = b</math>;
Rozwiąż <math>\displaystyle Ux = y</math>;
</pre>
 
Przypuśćmy, że taki rozkład <math>\displaystyle A=LU</math> istnieje. Wówczas, zapisując macierze w postaci
blokowej, eskponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy,
mamy
 
<center><math>\displaystyle
\begin{pmatrix}
a_{11} & a_{12}^T\\
a_{21} & A_{22}
\end{pmatrix}
=
\begin{pmatrix}
1 & 0^T\\
l_{21} & L_{22}
\end{pmatrix}
\begin{pmatrix}
u_{11} & u_{12}^T\\
0 & U_{22},
\end{pmatrix}
</math></center>
 
skąd (mnożąc blokowo macierz <math>\displaystyle L</math> przez <math>\displaystyle U</math>) wynika, że
* <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest
kopią pierwszego wiersza <math>\displaystyle A</math>,
* <math>\displaystyle l_{21} = a_{21}/u_{11}</math>, więc pierwsza kolumna <math>\displaystyle L</math> powstaje przez
podzielenie wszytkich elementów wektora <math>\displaystyle a_{21}</math> przez element na diagonali
<math>\displaystyle a_{11}</math>,
* <math>\displaystyle A_{22} - l_{21}u_{12}^T = L_{22}U_{22}</math>, a więc znalezienie podmacierzy
<math>\displaystyle L_{22}</math> oraz <math>\displaystyle U_{22}</math> sprowadza się do znalezienia rozkładu LU zmodyfikowanego
bloku <math>\displaystyle A_{22}</math> macierzy <math>\displaystyle A</math>,
wymiaru <math>\displaystyle (N-1)\times (N-1)</math>.
Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie
rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
klasycznych pętli. Jest to ważne w praktyce numerycznej, gdyż rekurencja
kosztuje: zarówno pamięć, jak i czas.
 
Ponadto zauważmy, że opisany algorytm możemy wykonać w miejscu, nadpisując
elementy <math>\displaystyle A</math> elementami macierzy <math>\displaystyle U</math> i <math>\displaystyle L</math> (jedynek z diagonali <math>\displaystyle L</math> nie musimy
pamiętać, bo wiemy a priori, że tam są).
 
<pre> [Algorytm rozkładu LU]
for k<nowiki>=</nowiki>1:N-1
if <math>\displaystyle a_{kk}</math> <nowiki>=</nowiki><nowiki>=</nowiki> 0
STOP;
end
for i<nowiki>=</nowiki>k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
<math>\displaystyle a_{ik}</math> <nowiki>=</nowiki> <math>\displaystyle a_{ik}/a_{ii}</math>;
end
for j<nowiki>=</nowiki>k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
for i<nowiki>=</nowiki>k+1:N
<math>\displaystyle a_{ij}</math> -<nowiki>=</nowiki> <math>\displaystyle a_{ik}a_{kj}</math>;
end
end
end
</pre>
 
Łatwo przekonać się, że <math>\displaystyle k</math>-ty obrót zewnętrznej pętli (tzn. <math>\displaystyle k</math>-ty krok
algorytmu rozkładu LU) kosztuje rzędu <math>\displaystyle 2(N-k)^2</math> operacji arytmetycznych, skąd łączny koszt tego
algorytmu rozkładu LU wynosi około <math>\displaystyle \frac{4}{3}N^3</math>.
 
Jeśli więc wykorzystać rozkład LU do rozwiązywania układu równań <math>\displaystyle Ax=b</math>, to mamy
następujące zestawienie kosztów:
* Koszt znalezienia rozkładu <math>\displaystyle A=LU</math>: <math>\displaystyle O(N^3)</math>;
* Koszt rozwiązania układu <math>\displaystyle Ly=b</math>: <math>\displaystyle O(N^2)</math>;
* Koszt rozwiązania układu <math>\displaystyle Ux=y</math>: <math>\displaystyle O(N^2)</math>.
Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania
wynosi już tylko <math>\displaystyle O(N^2)</math>.
 
{{uwaga|Złożoność obliczeniowa zadania rozwiązania układu równań liniowych||
 
Z powyższego wynika, że łączny koszt rozwiązania równania liniowego wynosi
<math>\displaystyle O(N^3)</math>. Można zastanawiać się, jaka jest ''najmniejsza możliwa'' liczba
operacji zmiennoprzecinkowych potrzebnych do rozwiązania układu równań
liniowych.
 
Można pokazać {Cormen}, że minimalny koszt rozwiązania układu <math>\displaystyle N</math> równań
liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch
macierzy <math>\displaystyle N\times N</math>. Znany jest całkiem prosty algorytm rekurencyjny,
wyznaczający iloczyn dwóch macierzy kosztem <math>\displaystyle 4.7\cdot N^{log_27} \approx 4.7
\cdot N^{2.807}</math> (algorytm Strassena). Bardziej skomplikowany (i praktycznie
nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt
<math>\displaystyle O(N^{2.376})</math>.  Tak więc równania liniowe daje się (teoretycznie) rozwiązać
kosztem <math>\displaystyle O(N^{2.376})</math>.
 
Jednak w praktyce nawet prosty algorytm Strassena zazwyczaj nie jest stosowany.
Wynika to stąd, że ma trochę gorsze własności numeryczne oraz wymaga dużo
dodatkowej pamięci na przechowywanie pośrednich wyników.
}}
 
===Wybór elementu głównego===
 
Opisany powyżej algorytm rozkładu LU niestety czasem może się załamać: gdy
napotka w czasie działania zerowy element na diagonali zmodyfikowanej
podmacierzy, np. chociaż macierz
 
<center><math>\displaystyle
A = \begin{pmatrix}  0 & 1\\ 1 & 0
\end{pmatrix}
</math></center>
 
jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od
razu zetknie się z dzieleniem przez <math>\displaystyle a_{11}=0</math>. Ale wystarczy zamienić ze sobą
kolejnością wiersze macierzy <math>\displaystyle A</math> (to znaczy, w układzie równań, zamienić ze sobą
miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez
problemu.
 
W praktyce obliczeniowej, aby uzyskać algorytm o [[<nowiki>|</nowiki>Uzupełnij]]{możliwie dobrych własnościach
numerycznych}, wykorzystujemy tzw. strategię ''wyboru elementu głównego w
kolumnie''.
Polega to na tym, że zanim wykonamy <math>\displaystyle k</math>-ty krok algorytmu rozkładu LU,
* szukamy w pierwszej kolumnie podmacierzy <math>\displaystyle A(k:N,k:N)</math> szukamy elementu o
największym module (taki element, na mocy założenia nieosobliwości macierzy,
jest niezerowy) --- to jest właśnie element główny
* zamieniamy ze sobą wiersz pierwszy  <math>\displaystyle A(k:N,k:N)</math> z wierszem, w którym
znajduje się element główny
* zapamiętujemy dokonaną permutację, bo potem --- gdy już przyjdzie do
rozwiązywania układu równań --- będziemy musieli dokonać
analogicznej permutacji wektora prawej strony
Wynikiem takiego algorytmu jest więc rozkład
 
<center><math>\displaystyle PA = LU,
</math></center>
 
gdzie <math>\displaystyle P</math> jest pewną (zerojedynkową) macierzą permutacji (tzn. macierzą
identyczności z przepermutowanymi wierszami).
 
Oprócz wyboru elementu głównego w kolumnie, stosuje się czasem inne strategie,
m.in. wybór w kolumnie oraz tzw. ''wybór pełny'', gdy elementu głównego szukamy w
''całej'' podmacierzy <math>\displaystyle A(k:N,k:N)</math>, co znacznie zwiększa liczbę porównań
niezbędnych do wskazania elementu głównego, ale też trochę poprawia własności
numeryczne takiego algorytmu.
 
W praktyce całą informację o wykonanych permutacjach przechowujemy nie w
zerojedynkowej macierzy j.w. ale w jednym wektorze.
 
<pre> [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
</pre>
 
Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
już rozkład LU wykonany z wyborem w kolumnie.
 
<pre> [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>;
</pre>
 
{}{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
 
<blockquote>  W optymalizacji szybkości działania programu numerycznego,
obecnie cała walka idzie o to, by procesor przez cały czas ''miał co liczyć''.
</blockquote>
 
Szczególnie jest to widoczne w algorytmach, które wykonują bardzo dużo operacji
na dużej liczbie
danych --- a tak jest m.in. w algorytmach algebry liniowej takich jak mnożenie
dwóch macierzy, czy rozwiązywanie układów równań liniowych: te algorytmy
najczęściej operują na <math>\displaystyle O(N^2)</math> danych i wykonują aż <math>\displaystyle O(N^3)</math> działań.
 
Ponieważ szybka pamięć komputerowa jest jednocześnie bardzo droga,
konstruktorzy komputerów osobistych stoją przed wykluczającymi się celami: z
jednej strony powinni zapewnić użytkownikowi jak najwięcej pamięci, z drugiej
zaś -- chcieliby zapewnić użytkownikowi jak najszybszą pamięć. Sprzeczność tych
dążeń ujawnia się, gdy dołączyć do nich trzecie, ale zasadnicze wymaganie:
całość ma być w rozsądnej cenie... Ostatecznie więc, z biegiem lat dramatycznie
pogłębia się przepaść pomiędzy  prędkością (podwajającą się, zgodnie z
heurystycznym prawem Moore'a, co półtora roku) procesora, a prędkością pamięci
RAM, do której procesor musi się odwoływać.
 
Powszechnie przyjętym sposobem pogodzenia tych sprzeczności jest pamięć
hierarchiczna.  Koncept polega na tym, że część pamięci, która najczęściej komunikuje się z  procesorem,
jest bardzo szybka (lecz jest jej relatywnie mało), natomiast pozostała pamięć
jest wolniejsza, za to może jej być bardzo dużo.
 
W praktycznej realizacji, zamiast dwóch poziomów pamięci (mała-szybka i
duża-wolna), w komputerze występuje wiele poziomów:
* rejestry procesora
* ''cache ''(pamięć podręczna) procesora
* ''cache ''drugiego poziomu (ostatnio także wbudowywana do procesora)
* pamięć operacyjna (RAM): główna pamięć komputera
* pamięć zewnętrzna (np. twardy dysk)
* pamięć masowa (CD-ROM, taśma streamera, itp.)
Efektywność działania komputera, zwłaszcza w wielkoskalowych obliczeniach
numerycznych, bardzo istotnie zależy od skutecznego wykorzystania hierarchii
pamięci tak, by jak największa część operacji była wykonywana na zmiennych
znajdujących się w danej chwili w jak najszybszej pamięci procesora.
 
Kluczem do skutecznego wykorzystania hierarchii pamięci jest zasada
lokalności w czasie i w przestrzeni:
 
<blockquote> 
; Lokalność w czasie:
:  Używać danego fragmentu pamięci intensywnie, ale rzadko.
 
; Lokalność w przestrzeni (adresowej):
:  W danej chwili, odnosić
się do adresów pamięci leżących blisko siebie.
</blockquote>
 
Zachowanie zasady lokalności w czasie jest bardzo ważne dla efektywnego
wykorzystania cache'a -- ze względu na małe rozmiary tej pamięci, wymiana
zmiennych może być tam intensywna. Zasada lokalności w przestrzeni też jest
ważna, przy czym nie tylko ze względu na efektywne wykorzystanie cache'a, ale
także dla efektywnego wykorzystania pamięci
wirtualnej.
 
===Jak napisać kod źle 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 <span style="font-style:monospace"> double </span> 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====
 
<pre> [ijk]
/* ijk */
for (i <nowiki>=</nowiki> 0; i < N; i++)
for (j <nowiki>=</nowiki> 0; j < N; j++)
for (k <nowiki>=</nowiki> 0; k < N; k++)
C[i*N+j] +<nowiki>=</nowiki> A[i*N+k]*B[k*N+j];
</pre>
 
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 <span style="font-style:monospace"> double </span>.
 
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 <span style="font-style:monospace"> B[k*N+j] </span>, <span style="font-style:monospace"> k </span> <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:
 
<pre> [ikj]
/* ikj */
for (i <nowiki>=</nowiki> 0; i < N; i++)
for (k <nowiki>=</nowiki> 0; k < N; k++)
for (j <nowiki>=</nowiki> 0; j < N; j++)
C[i*N+j] +<nowiki>=</nowiki> A[i*N+k]*B[k*N+j];
</pre>
 
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 <span style="font-style:monospace"> B </span>, 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  <span style="font-style:monospace"> B </span>. 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>:
 
<pre> [bikj(16)]
/* bikj(16) */
for (i <nowiki>=</nowiki> 0; i < N; i+<nowiki>=</nowiki>16)
for (k <nowiki>=</nowiki> 0; k < N; k+<nowiki>=</nowiki>16)
for (j <nowiki>=</nowiki> 0; j < N; j+<nowiki>=</nowiki>16)
for (ii <nowiki>=</nowiki> i; ii < i+15; ii++)
for (kk <nowiki>=</nowiki> k; kk < k+15; kk++)
for (jj <nowiki>=</nowiki> j; jj < j+15; jj++)
C[ii*N+jj] +<nowiki>=</nowiki> A[ii*N+kk]*B[kk*N+jj];
 
</pre>
 
(podobnie działał algorytm bikj(32), tym razem na blokach <math>\displaystyle 32\times 32</math>).
 
Wadą poprzedniego algorytmu, ikj, było to, że w dwu zewnętrznych pętlach powracał do
wszystkich <math>\displaystyle N^2</math> wartości <math>\displaystyle C</math> i <math>\displaystyle B</math>, przecząc zasadzie ''lokalności  w
czasie''. Wydzielenie w algorytmie bijk(16) operacji na blokach macierzy ma na celu uniknięcie tego
problemu: trzy najbardziej zewnętrzne pętle to <math>\displaystyle (1024/16)^3</math> obrotów, czyli
czterokrotnie mniej niż dwie najbardziej zewnętrzne pętle w poprzednim
algorytmie. Gdyby udało się zachować liczbę ''cache misses ''na poprzednim poziomie,
można byłoby liczyć na czterokrotne przyspieszenie. ('''Do sprawdzenia: I tak z grubsza jest:
teoretycznie wszystkie 3 podmacierze powinny mieścić się w ''cache '''u.''')
 
====Algorytmy DGEMM i ATLAS DGEMM====
 
Algorytm DGEMM to był algorytm mnożenia macierzy z pakietu BLAS --- jest to
właśnie profesjonalny algorytm blokowy, ale niezoptymalizowany na naszą
architekturę. Dopiero algorytm DGEMM podrasowany w pakiecie ATLAS dał nam
sprawność wynoszącą 1.2 Gflopów na sekundę -- czyli praktycznie
maksimum (''teoretycznie'', z Durona można wycisnąć dwa flopy w cyklu
zegara, co dawałoby <math>\displaystyle r_{\max}</math> <nowiki>=</nowiki> 2.2 Gflop/s, ale w praktyce jest to mało
prawdopodobne) tego,
co da się wycisnąć z tej wersji Durona na liczbach podwójnej precyzji.
 
===Macierze w pamięci komputera===
 
Do implementacji algorytmów numerycznych powszechnie używa się dwóch języków:
Fortranu i C. Zgodnie z naszą konwencją, skupimy się na programach w C, nie
możemy jednak przejść obojętnie wobec faktu, że jest bardzo wiele znakomitego
oprogramowania numerycznego w Fortranie. W Rozdziale&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>
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<nowiki>|</nowiki>Uzupełnij]]{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 <span style="font-style:monospace"> matrix </span>, gdyż pętle są
ustawione tak, by przechodzić przez macierz wzdłuż kolejnych kolumn. Dlatego nie
jest tu konieczne użycie makra <span style="font-style:monospace"> IJ() </span>, a sprytne wykorzystanie
pointera <span style="font-style:monospace"> ptr </span>
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 <span style="font-style:monospace"> IJ() </span> (i obliczania wyrażenia <span style="font-style:monospace"> i+j*N </span>) 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 <span style="font-style:monospace"> IJ() </span>  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 <span style="font-style:monospace"> j*N </span>, 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 <span style="font-style:monospace"> IJ() </span>) 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|]].
 
===Włączenie zewnętrznej biblioteki fortranowskiej
do programu===
 
Wiele spośród doskonałych bibliotek numerycznych zostało napisanych w Fortranie
77 (np. ARPACK, LAPACK, ODEPACK). Tymczasem, nasze programy zdecydowaliśmy się
(ze względów wymienionych w Rozdziale&nbsp;[[##sec:|Uzupelnic sec:|]]) pisać w języku C. Na
szczęście, istnieje prosty sposób wykorzystania gotowych pakietów
fortranowskich przez zewnętrzne programy w C: wystarczy skorzystać z genialnej
biblioteki <span style="font-style:monospace"> f2c </span> lub jej modyfikacji na użytek kompilatora GCC,
biblioteki <span style="font-style:monospace"> gfortran </span>.
 
Najczęściej jest tak, że daną bibliotekę (fortranowską) instalujemy na swoim
komputerze z plików źródłowych (np. ściągniętych z Internetu). Instalacja
takiej biblioteki, powiedzmy, LAPACK'a, kończy się utworzeniem pliku
<span style="font-style:monospace"> liblapack.a </span>, zawierającego skompilowane wszystkie funkcje tej
biblioteki.
 
Z uwagi na względnie powszechną dostępność LAPACKa w pakietach RPM,
właśnie na przykładzie tych bibliotek omówimy sposób wykorzystania bibliotek
fortranowskich w
programie w C.
 
Napiszemy program, który będzie liczył normę zadanego
wektora, korzystając z funkcji <span style="font-style:monospace"> DNRM2 </span> biblioteki BLAS.
 
Najpierw musimy zorientować się, jak wygląda schemat wywołania takiej funkcji w
Fortranie. Otóż funkcja wygląda następująco:
 
      DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX )
*    .. Scalar Arguments ..
      INTEGER                          INCX, N
*    .. Array Arguments ..
      DOUBLE PRECISION                  X( * )
*    ..
*
*  DNRM2 returns the euclidean norm of a vector via the function
*  name, so that
*
*    DNRM2 :<nowiki>=</nowiki> sqrt( x'*x )
*
i tak dalej...
Tak więc nasza funkcja obliczająca normę wektora ma trzy argumenty: <span style="font-style:monospace"> N </span> --
długość wektora (<span style="font-style:monospace"> INTEGER </span>), <span style="font-style:monospace"> X </span> -- wektor, którego długość chcemy
obliczyć (tablica liczb <span style="font-style:monospace"> DOUBLE PRECISION </span>) oraz tajemniczy dodatkowy parametr
<span style="font-style:monospace"> INCX </span> typu <span style="font-style:monospace"> INTEGER </span> -- jest to wartość skoku, określająca co który
element wektora uwzględnić w obliczaniu normy: aby policzyć normę całego
wektora, bierzemy <span style="font-style:monospace"> INCX </span> równe 1. Używając zapisu Octave, <span style="font-style:monospace"> DNRM2 </span>
oblicza po prostu
 
<pre>
norm( X(1:INCX:N) )
</pre>
 
Kod obiektowy tej funkcji znajduje się już w bibliotece BLAS, zawartej w pliku
<span style="font-style:monospace"> libblas.a </span>. Chcielibyśmy wykorzystać tę funkcję w programie w C, a jak
wiadomo, każda funkcja w C powinna mieć swój prototyp. Jaki powinien być ''prototyp'' tej funkcji?
 
Przede wszystkim, zacznijmy od nazwy. W przypadku kompilatora
<span style="font-style:monospace"> gcc </span>/<span style="font-style:monospace"> gfortran </span>, nazwą funkcji do wykorzystania w C będzie
<span style="font-style:monospace"> dnrm2_ </span> (tak! małymi literami i z przyrostkiem "<span style="font-style:monospace"> _ </span>").
 
Jeśli zaś chodzi o argumenty, to zapewne co do drugiego nie będziemy mieli
wątpliwości: jako wektor <span style="font-style:monospace"> X </span> przekażemy -- naturalnie -- ''wskaźnik'' do
tablicy <span style="font-style:monospace"> X </span> (typu <span style="font-style:monospace"> double </span>), czyli po prostu: jej nazwę. Co z
pozostałymi argumentami? Okazuje się, że reguła jest niezmienna:
 
<blockquote>  Każdy argument funkcji fortranowskiej zastępujemy ''wskaźnikiem'' do odpowiedniego typu:
 
{| border=1
|+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
|-
|
Fortran 77  ||  C
|-
|
INTEGER  ||  int
|-
| REAL  ||  float
|-
| DOUBLE PRECISION  ||  double
|-
| COMPLEX  ||  struct { float Re, Im; }
|-
| DOUBLE COMPLEX  ||  struct { double Re, Im; }
|-
| CHARACTER  ||  char
|-
|
 
|}
 
</blockquote>
 
<blockquote>  Wszystkim argumentom macierzowym danego typu w Fortranie
(reprezentującym macierze jedno-, dwu-, i więcejwymiarowe) przypisujemy w C
(pojedynczy) wskaźnik do tego typu (o czym w będzie mowa w następnym
przykładzie). </blockquote> 
 
A więc pierwszym i trzecim argumentem funkcji <span style="font-style:monospace"> dnrm2_ </span>
będą wskaźniki do <span style="font-style:monospace"> int </span>. Ponieważ
funkcja <span style="font-style:monospace"> DNRM2 </span> zwraca w wyniku liczbę podwójnej precyzji, to ostatecznie
prototyp naszej funkcji w C będzie następujący:
 
<pre>
double dnrm2_(int* N, double* X, int* INCX);
</pre>
 
No to wykorzystajmy naszą funkcję:
 
<pre> [Wykorzystanie funkcji DNRM2 fortranowskiej biblioteki BLAS w programie w C]
#include <stdio.h>
double dnrm2_(int*,double*,int*);
 
int main(void)
{
int n, incx<nowiki>=</nowiki>1;
double x[3]<nowiki>=</nowiki> {0,1,2};
n <nowiki>=</nowiki> 3;
printf("Norma podanego wektora:
return(0);
}
</pre>
 
Zwróćmy uwagę na sposób kompilacji tego programu:
 
gcc -o testblas testblas.c -lblas -lgfortran -lm
oprócz biblioteki BLAS, co
naturalne, musimy dołączyć jeszcze bibliotekę matematyczną (może być
wykorzystywana przez BLAS!) oraz, co bardzo ważne,  specjalną bibliotekę:
<span style="font-style:monospace"> gfortran </span>, umożliwiającą koegzystencję Fortranu i
C.
 
====Funkcja fortranowska z argumentem macierzowym====
 
Należy szczególną uwagę zwrócić na argumenty macierzowe funkcji w Fortranie,
gdyż bardzo łatwo tu o pomyłkę, którą potem trudno wychwycić. Aby niepotrzebnie
nie komplikować przykładu subtelnościami funkcji BLAS, rozważmy kod źródłowy
prostej funkcji w Fortranie 77, która po prostu wypełnia liczbami pewną macierz
wymiaru <math>\displaystyle M\times N</math>:
 
[Funkcja w Fortranie 77, wypełniająca macierz <math>\displaystyle M\times N</math>]
SUBROUTINE FILLMATRIX ( M, N, MATRIX )
INTEGER M,N
DOUBLE PRECISION MATRIX(M, N)
DO 10 I<nowiki>=</nowiki>1,M
DO 20 J<nowiki>=</nowiki>1,N
MATRIX(I,J) <nowiki>=</nowiki> I+2*J
20 CONTINUE
10 CONTINUE
END
Nawet osoby nie znające Fortranu domyślą się, że wynikiem działania naszej
funkcji, np. dla <math>\displaystyle M=2</math>, <math>\displaystyle N=5</math>, będzie macierz
 
<center><math>\displaystyle
\lstF{MATRIX} =
\begin{pmatrix}
3 & 5 & 7 & 9 & 11 \\
4 & 6 & 8 & 10 & 12
\end{pmatrix}
</math></center>
 
Naturalnie, możemy wywołać ją wprost z programu w C przy użyciu poprzednio
poznanej techniki i następującego kodu (tym razem prototyp funkcji
<span style="font-style:monospace"> fillmatrix_ </span> umieszczamy w osobnym pliku nagłówkowym <span style="font-style:monospace"> ffortran.h </span>,
gdzie mamy zamiar kolekcjonować prototypy w C lokalnie wykorzystywanych funkcji
fortranowskich):
 
<pre> [Wykorzystanie funkcji fortranowskiej
operującej na macierzy. Zwróćmy uwagę na sposób użycia argumentu
macierzowego]
#include <stdio.h>
#include <stdlib.h>
int fillmatrix_(int *, int *, double *);
 
int main()
{
int MM, NN, i, j;
double *A;
  MM <nowiki>=</nowiki> 2; NN <nowiki>=</nowiki> 5;
A <nowiki>=</nowiki> (double *)malloc(MM*NN*sizeof(double));
 
  fillmatrix_( &MM, &NN, A );
 
printf(" elementy wektora A:"); 
  for ( i <nowiki>=</nowiki> 0; i < NN*MM ; i++ ){
printf("}
 
printf(" A zinterpretowany jako macierz:"); 
for ( j <nowiki>=</nowiki> 0 ; j < MM ; j++ )
{
  for ( i <nowiki>=</nowiki> 0; i < NN ; i++ )
  printf("printf("");
  }
 
  free( A );
 
  return(0);
}
</pre>
 
Zauważmy, że mimo iż funkcja fortranowska odwołuje się do macierzy ''dwuwymiarowej'', jako argument jej wywołania w C przekazujemy ''tablicę
jednowymiarową'' odpowiedniej wielkości.
 
===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 <span style="font-style:monospace"> float </span> i <span style="font-style:monospace"> double </span> i około pięciokrotne  na zmiennych
typu <span style="font-style:monospace"> complex  </span> i <span style="font-style:monospace"> double complex </span>.
 
Aby osiągnąć największe przyspieszenie, bibliotekę ATLAS należy skompilować
samemu na własnej (''nieobciążonej'' w trakcie kompilacji!) architekturze. W plikach
<span style="font-style:monospace"> Makefile </span> ATLASa brak jednak opcji instalacji bibliotek w standardowych
miejscach --- trzeba zrobić to samemu.
 
BLAS i LAPACK są często wykorzystywane przez inne biblioteki numeryczne,
na nich opierają się również funkcje macierzowe w Octave i MATLABie.
Optymalizowane biblioteki BLAS i LAPACK dla swoich architektur oferują
producenci procesorów Intel (biblioteka MKL) oraz AMD (biblioteka ACML)
 
BLAS {BLAS-home-page} jest kolekcją procedur służących manipulacji podstawowymi obiektami
algebry liniowej: skalarami, wektorami i macierzami. Obsługiwane są obiekty
rzeczywiste (w pojedynczej albo podwójnej precyzji) oraz zespolone (także w
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;
* 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
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ń
arytmetycznych przy
<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.
Wówczas dany algorytm algebry liniowej najlepiej zapisać przy użyciu procedur  BLAS Level
3, naturalnie, pod warunkiem, że w ogóle daje się to zrobić; typową strategią w
takim wypadku jest tworzenie algorytmów blokowych, operujących na ''blokach''
macierzy, zamiast na jej pojedynczych elementach.
 
Większość użytkowników BLAS nie będzie jednak miała potrzeby pisania własnych
algorytmów blokowych, gdyż funkcje rozwiązujące podstawowe zadania numerycznej
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 <span style="font-style:monospace"> PRRFF </span>, gdzie
;
:  <span style="font-style:monospace"> P </span> 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;
 
;
:  <span style="font-style:monospace"> RR </span> oznacza rodzaj zadania, np. GE oznacza ''GEneral '', czyli zadanie ogólne
(praktycznie bez założeń), a SY oznacza ''SYmmetric '', czyli zadanie symetryczne;
 
;
:  <span style="font-style:monospace"> FF </span> 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.
 
Naszym zadaniem jest wykonanie operacji
<center><math>\displaystyle
y \leftarrow  \alpha A x + y,
</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>
współrzędnych.
 
To zadanie realizuje procedura BLASów o nazwie
<span style="font-style:monospace"> DGEMV </span>. W rzeczywistości, ta procedura wykonuje ogólniejsze zadanie
wyznaczania wektora
 
<center><math>\displaystyle
y \leftarrow  \alpha B x + \beta y,
</math></center>
 
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
razem argumentem macierzowym, jaki przekazujemy <span style="font-style:monospace"> DGEMV </span>, jest wyjściowa
macierz <math>\displaystyle A</math>).
 
Jak wiemy, że jest możliwe bezpośrednie wykorzystanie
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 <span style="font-style:monospace"> DGEMV </span> ze
strony {opis:dgemv}, w programie w C powinniśmy
napisać prototyp tej funkcji następująco:
 
<pre>
int dgemv_( char* TRANS,
int* M,
int* N,
double* ALPHA,
double* A, 
int* LDA, 
double* X, 
int* INCX,
double* BETA,
double* Y,
int* INCY );
</pre> 
 
Dla własnej wygody, a także dla przyszłego wykorzystania, umieścimy ten
prototyp, razem z innymi przydatnymi drobiazgami (m.in. makro <span style="font-style:monospace"> IJ </span>
dające wygodny dostęp do macierzy niezależny od jej wewnętrznej reprezentacji, a
także zmienne całkowite
<span style="font-style:monospace"> static int BLASONE <nowiki>=</nowiki> 1, BLASMONE <nowiki>=</nowiki> -1; </span>), w pliku
nagłówkowym <span style="font-style:monospace"> blaslapack.h </span>.
 
Dalej już jest łatwo: oto pełny kod programu realizującego operację mnożenia macierzy przez wektor
przy użyciu procedury BLAS <span style="font-style:monospace"> DGEMV </span>:
 
<pre>
#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);
}
</pre>
 
Zwróćmy uwagę na wykorzystanie "stałej" <span style="font-style:monospace"> BLASONE </span>, równej 1,
predefiniowanej w pliku <span style="font-style:monospace"> blaslapack.h </span>. 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 <span style="font-style:monospace"> -llapack </span> 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.
Dla osiągnięcia maksymalnej efektywności kodu, trzeba skorzystać z
optymalizowanych BLAS, które obecnie są dostępne nawet w kilku wariantach na
architektury x86.

Wersja z 13:17, 28 sie 2006

Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki

Algorytm

UULLLULA 

<math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
for (i = N-1; i >= 1; i--)
	<math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
	u_{i,j}x_j^*\right)/u_{i,i}</math>;
Źródło: „https://wazniak.mimuw.edu.pl/index.php?title=MN04&oldid=31492”