MN04

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

$Ax=b,$

gdzie $A$ jest nieosobliwą macierzą $N\times N$, a dany wektor prawej strony $b\in R^N$.

W praktyce spotyka się zadania z $N=2,3,\dots 1000$. Zdarzają się także czaem specjalne macierze o wymiarach nawet rzędu $10^8$!

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 $A^{-1}$ i następnie $x=A^{-1}b$

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ą

Rozważmy układ z macierzą trójkątną $A$. Będą nas szczególnie interesować macierze trójkątne górne, dla których $a_{i,j}=0$ gdy $i>j$, oraz macierze trójkątne dolne z jedynkami na przekątnej, tzn. $a_{i,j}=0$, $i<j$, oraz $a_{i,i}=1$. Macierze pierwszego rodzaju będziemy oznaczać przez $U$, a drugiego rodzaju przez $L$.

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \displaystyle L = \beginpmatrix 1 & & & & & \\ * & 1 & & & & \\ * & * & 1 & & & \\ * & * & * & 1 & & \\ \vdots & \vdots & \vdots & \ddots & \ddots & \\ * & * & * & \cdots & * & 1 \endpmatrix , \qquad U = \beginpmatrix * & * & * & * & \cdots & * \\ & * & * & * & \cdots & * \\ & & * & * & \cdots & * \\ & & & * & \ddots & \vdots \\ & & & & \ddots & * \\ & & & & & * \endpmatrix }

Układ z nieosobliwą macierzą trójkątną górną

$U\overrightarrow{x}=\overrightarrow{c},$

$U=(u_{i,j})$, $\overrightarrow{c}=(c_j)$, można rozwiązać stosując algorytm:

Algorytm

---

$x_N^{*}=c_N/u_{N,N}$; for (i = N-1; i >= 1; i--) $x_i^{*}:=(c_i-\sum _{j=i+1}^N{u}_{i,j}{x}_j^{*})/u_{i,i}$;

(Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy implikuje, że $u_{i,i}\ne 0$, $\mathrm{\forall }i$.) Podobnie, układ $L\overrightarrow{x}=\overrightarrow{c}$ rozwiązujemy algorytmem:

Algorytm

---

$x_1=c_1$; for (i=2; i <= N; i++) $x_i=c_i-\sum _{j=1}^{i-1}{l}_{i,j}{x}_j^{*}$;

Oba algorytmy wymagają rzędu $N^2/2$ mnożeń lub dzieleń i $N^2/2$ dodawań lub odejmowań, a więc łącznie $O(N^2)$ działań arytmetycznych.

Układy z macierzą ortogonalną

Równie tanio można rozwiązać układ równań

$Q\overrightarrow{x}=\overrightarrow{b},$

gdy $Q$ jest macierzą ortogonalną, to znaczy $Q^{T}Q=I$. Rzeczywiście, z ortogonalności mamy natychmiast, że

$\overrightarrow{x}=Q^T\overrightarrow{b}$

i w konsekwencji $x$ można wyznaczyć kosztem takim jak koszt mnożenia macierzy przez wektor, czyli $O(N^2)$ operacji.

Podobnie, gdy $Q\in C^{N\times N}$ jest unitarna, to znaczy $Q^{H}Q=I$, rozwiązaniem układu równań jest

$\overrightarrow{x}=Q^H\overrightarrow{b}.$

Metoda eliminacji Gaussa

W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym rozwiązywania układu równań

$Ax=b$

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 $L$ (z jedynkami na diagonali) oraz trójkątnej górnej $U$ takich, że

$A=LU,$

a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:

Algorytm Algorytm eliminacji Gaussa

---

 AEG
Znajdź rozkład $A=LU$;
Rozwiąż $Ly=b$;
Rozwiąż $Ux=y$;

Przypuśćmy, że taki rozkład $A=LU$ istnieje. Wówczas, zapisując macierze w postaci blokowej, eskponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy, mamy

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \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 }

skąd (mnożąc blokowo macierz $L$ przez $U$) wynika, że

• $u_{11}=a_{11}$ oraz $u_{12}=a_{12}$, więc pierwszy wiersz $U$ jest

kopią pierwszego wiersza $A$,

• $l_{21}=a_{21}/u_{11}$, więc pierwsza kolumna $L$ powstaje przez

podzielenie wszytkich elementów wektora $a_{21}$ przez element na diagonali $a_{11}$,

• $A_{22}-l_{21}{u}_{12}^T=L_{22}{U}_{22}$, a więc znalezienie podmacierzy

$L_{22}$ oraz $U_{22}$ sprowadza się do znalezienia rozkładu LU zmodyfikowanego bloku $A_{22}$ macierzy $A$, wymiaru $(N-1)\times (N-1)$.

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 $A$ elementami macierzy $U$ i $L$ (jedynek z diagonali $L$ nie musimy pamiętać, bo wiemy a priori, że tam są).

Algorytm title=Algorytm rozkładu LU

---

for k=1:N-1 if $a_{kk}$ == 0 STOP; end for i=k+1:N /* wyznaczenie $k$-tej kolumny $L$ */ $a_{ik}$ = $a_{ik}/a_{ii}$; end for j=k+1:N /* modyfikacja podmacierzy $A_{22}=A_{22}-l_{21}{u}_{12}^T$ */ for i=k+1:N $a_{ij}$ -= $a_{ik}{a}_{kj}$; end end end

Łatwo przekonać się, że $k$-ty obrót zewnętrznej pętli (tzn. $k$-ty krok algorytmu rozkładu LU) kosztuje rzędu $2(N-k{)}^2$ operacji arytmetycznych, skąd łączny koszt tego algorytmu rozkładu LU wynosi około $\frac{4}{3}{N}^3$.

Jeśli więc wykorzystać rozkład LU do rozwiązywania układu równań $Ax=b$, to mamy następujące zestawienie kosztów:

1. Koszt znalezienia rozkładu $A=LU$: $O(N^3)$;
2. Koszt rozwiązania układu $Ly=b$: $O(N^2)$;
3. Koszt rozwiązania układu $Ux=y$: $O(N^2)$.

Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania wynosi już tylko $O(N^2)$.

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 $O(N^3)$. 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 $N$ równań liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch macierzy $N\times N$. Znany jest całkiem prosty algorytm rekurencyjny, wyznaczający iloczyn dwóch macierzy kosztem $4.7\cdot N^{lo{g}_{2}7}\approx 4.7\cdot N^{2.807}$ (algorytm Strassena). Bardziej skomplikowany (i praktycznie nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt $O(N^{2.376})$. Tak więc równania liniowe daje się (teoretycznie) rozwiązać kosztem $O(N^{2.376})$.

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

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \displaystyle A = \beginpmatrix 0 & 1\\ 1 & 0 \endpmatrix }

jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od razu zetknie się z dzieleniem przez $a_{11}=0$. Ale wystarczy zamienić ze sobą kolejnością wiersze macierzy $A$ (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 $k$-ty krok algorytmu rozkładu LU,

1. szukamy w pierwszej kolumnie podmacierzy $A(k:N,k:N)$ 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

1. zamieniamy ze sobą wiersz pierwszy $A(k:N,k:N)$ z wierszem, w którym

znajduje się element główny

1. 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

$PA=LU,$

gdzie $P$ 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 $A(k:N,k:N)$, 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=Algorytm rozkładu LU z wyborem elementu głównego w kolumnie

---

P = 1:N; /* tu zapamiętujemy wykonane permutacje */ for k=1:N-1

w wektorze A(k:N,k) znajdź element główny $a_{pk}$; zamień ze sobą wiersze A(k,k:N) i A(p,k:N); P(k) = p; P(p) = k; if $a_{kk}$ STOP: macierz osobliwa! end

/* kontunuuj tak jak w algorytmie bez wyboru */

for i=k+1:N /* wyznaczenie $k$-tej kolumny $L$ */ $a_{ik}$ = $a_{ik}/a_{ii}$; end for j=k+1:N /* modyfikacja podmacierzy $A_{22}=A_{22}-l_{21}{u}_{12}^T$ */ for i=k+1:N $a_{ij}$ -= $a_{ik}{a}_{kj}$; 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=Algorytm eliminacji Gaussa z wyborem elementu głównego w kolumnie

---

znajdź rozkład $PA=LU$; rozwiąż względem $y$ układ z macierzą górną trójkątną $Uy=Pb$; rozwiąż względem $x$ układ z macierzą dolną trójkątną $Lx=y$;

{}{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: $A=A^T$ oraz $x^{T}Ax>0$, $\mathrm{\forall }x\ne 0$,
• silnie diagonalnie dominujących: macierz $A$ (lub $A^T$) spełnia

$|a_{ii}|>\sum _{j\ne i}|a_{ij}|,\mathrm{\forall }i.$

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 $O(N^2)$ danych i wykonują aż $O(N^3)$ 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 $1024\times 1024$ przy użyciu kilku matematycznie równoważnych algorytmów (nazwaliśmy je umownie ijk, ikj, bikj($\cdot$) --- 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ł 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ł 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:

Uzupelnij tytul

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 = 0; i < N; i++) for (j = 0; j < N; j++) for (k = 0; k < N; k++) C[i*N+j] += 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 $A$ oraz $B$ powodujemy, że przy odwoływaniu się do $B$, {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! = $0\dots N$, odwołujemy się do co 1024. elementu wektora B, a zatem kolejne odwołania lądują w tej samej sekcji {cache}'a -- która mieści ledwie 2 linie. Skutkiem tego, że w każdym obrocie pętli mamy chybienie, jest złudzenie pracowania z komputerem bez pamięci {cache} (a nawet gorzej, bo {cache miss} dodatkowo kosztuje) czyli jakby z komputerem o szybkości 10 MHz = 100 MHz/10 (bo magistrala ({bus}) jest taktowana 100 MHz, a odwołanie do pamięci RAM kosztuje mniej więcej 10 cykli magistrali). I rzeczywiście, wyniki zdają się to potwierdzać.

Algorytm ikj

Różni się on od poprzedniego jedynie kolejnością dwóch wewnętrznych pętli:

Algorytm

---

/* ikj */ for (i = 0; i < N; i++) for (k = 0; k < N; k++) for (j = 0; j < N; j++) C[i*N+j] += A[i*N+k]*B[k*N+j];

Okazuje się, że taka prosta zmiana dramatycznie poprawia sytuację!

Tym razem, w odwołaniu do $B$ 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 $16\times 16$:

Algorytm

---

/* bikj(16) */ for (i = 0; i < N; i+=16) for (k = 0; k < N; k+=16) for (j = 0; j < N; j+=16) for (ii = i; ii < i+15; ii++) for (kk = k; kk < k+15; kk++) for (jj = j; jj < j+15; jj++) C[ii*N+jj] += A[ii*N+kk]*B[kk*N+jj];

(podobnie działał algorytm bikj(32), tym razem na blokach $32\times 32$).

Wadą poprzedniego algorytmu, ikj, było to, że w dwu zewnętrznych pętlach powracał do wszystkich $N^2$ wartości $C$ i $B$, 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 $(1024/16{)}^3$ 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 $r_{\max }$ = 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 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ą $n\times m$ o elementach $a_{ij}$, $i=1\dots n$, $j=1\dots m$,

Parser nie mógł rozpoznać (nieznana funkcja „\beginpmatrix”): {\displaystyle \displaystyle \beginpmatrix a_{11} & \cdots & a_{1m}\\ \vdots & & \vdots\\ a_{n1} & \cdots & a_{nm} \endpmatrix . }

to kolejne miejsca w przestrzeni adresowej zajmują elementy

$a_{11},a_{21},\dots ,a_{n1},a_{12},a_{22},\dots ,a_{n2},\dots a_{nm}.$

Dla odmiany, C przechowuje w pamięci elementy macierzy wierszami, tzn. kolejno

$a_{11},a_{12},\dots ,a_{1m},a_{21},a_{22},\dots ,a_{2m},\dots a_{nm}.$

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:

1. przydzielić pamięć na tablicę wskaźników do wierszy
2. 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 $n\times m$ będziemy zapisywali do wektora o długości $n\cdot m$.

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 $(i,j)$-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 $a_{ij}$, $i=1\dots n$, $j=1\dots m$.

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ł 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ć 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, 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

$y\leftarrow \alpha x+y,$

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

$y\leftarrow \alpha Ax+y$

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ł Uzupelnic sec:architektura|;

• BLAS Level 3 -- operacje typu macierz--macierz, np. mnożenie dwóch

macierzy:

$C\leftarrow \alpha A\cdot B+C$

W związku z tym, że operacje macierzowe wymagają wykonania $O(N^3)$ działań arytmetycznych przy $O(N^2)$ danych (gdzie $N$ jest wymiarem macierzy), wykorzystanie zoptymalizowanych procedur BLAS Level 3 może znacząco przyspieszyć wykonywanie obliczeń na maszynach z pamięcią hierarchiczną (zob. Rozdział 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 {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.

Uzupelnij tytul

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$T^{}$ 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

Przykład 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 Uzupelnic sec:FortranC|.

Naszym zadaniem jest wykonanie operacji

$y\leftarrow \alpha Ax+y,$

gdzie $A$ jest zadaną macierzą $N\times M$, natomiast $y$ jest wektorem o $M$ współrzędnych.

To zadanie realizuje procedura BLASów o nazwie {DGEMV}. W rzeczywistości, ta procedura wykonuje ogólniejsze zadanie wyznaczania wektora

$y\leftarrow \alpha Bx+\beta y,$

przy czym macierz $B$ może być równa albo $A$, albo $A^T$ (jednak za każdym razem argumentem macierzowym, jaki przekazujemy {DGEMV}, jest wyjściowa macierz $A$).

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 {DGEMV} ze strony {opis:dgemv}, w programie w C powinniśmy napisać prototyp tej funkcji następująco:

Algorytm

---

int dgemv_( char* TRANS, int* M, int* N, 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ł Uzupelnic sec:macierze-w-komputerze| i zmienne całkowite !static int BLASONE = 1, BLASMONE = -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

---

#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>= mmread("mbeacxc.mtx", &N, &M );

x = (double *)malloc(N*sizeof(double)); y = (double *)malloc(M*sizeof(double)); for (i = 1; i <= N; i++) x[IJ(i,1,N)] = (double)i;

/* obliczamy y = 5*A*x, korzystając z procedury BLAS L2: DGEMV */

{ char TRANS = 'N'; double ALPHA = 5.0, BETA = 0.0;

dgemv_(&TRANS, &N, &M, &ALPHA, A, &N, x, &BLASONE, &BETA, y, &BLASONE );

}

/* wydruk wyniku */ for (i = 1; i <= 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!