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

$L=\left(\begin{array}{cccccc}1& & & & & \\ *& 1& & & & \\ *& *& 1& & & \\ *& *& *& 1& & \\ ⋮& ⋮& ⋮& \ddots & \ddots & \\ *& *& *& \cdots & *& 1\end{array}\right),U=\left(\begin{array}{cccccc}*& *& *& *& \cdots & *\\ & *& *& *& \cdots & *\\ & & *& *& \cdots & *\\ & & & *& \ddots & ⋮\\ & & & & \ddots & *\\ & & & & & *\end{array}\right)$

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

$Ux=c,$

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

Algorytm Podstawienie w tył

---

<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>;

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

Algorytm Podstawienie w przód

---

<math>\displaystyle x_1 = c_1</math>;
for (i=2; i <= N; i++)
	<math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;

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ń

$Qx=b,$

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

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

$x=Q^{H}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 Rozwiązanie układu z wykorzystaniem rozkładu LU

---

Znajdź rozkład <math>\displaystyle A=LU</math>;
Rozwiąż <math>\displaystyle Ly = b</math> przez podstawienie w przód;
Rozwiąż <math>\displaystyle Ux = y</math> przez podstawienie w tył;

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

$\left(\begin{array}{cc}{a}_{11}& a_{12}^T\\ a_{21}& A_{22}\end{array}\right)=\left(\begin{array}{cc}1& 0^T\\ l_{21}& L_{22}\end{array}\right)\left(\begin{array}{cc}{u}_{11}& u_{12}^T\\ 0& U_{22},\end{array}\right)$

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 Rozkład LU metodą eliminacji Gaussa

---

for k=1:N-1
	if <math>\displaystyle a_{kk}</math> == 0 
		STOP;
	end
	for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
		<math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
	end
	for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
		for i=k+1:N 
			<math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
		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:

• Koszt znalezienia rozkładu $A=LU$: $O(N^3)$;
• Koszt rozwiązania układu $Ly=b$: $O(N^2)$;
• 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ć \cite{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

$A=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$

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 [[|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 $k$-ty krok algorytmu rozkładu LU,

• 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

• zamieniamy ze sobą wiersz pierwszy $A(k:N,k:N)$ 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

$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 Rozkład 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 <math>\displaystyle a_{pk}</math>;
	zamień ze sobą wiersze A(k,k:N) i A(p,k:N);
	P(k) = p; P(p) = k;
	if <math>\displaystyle a_{kk}</math>
		STOP: macierz osobliwa!
	end
	
	/* kontunuuj tak jak w algorytmie bez wyboru */
	
	for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
		<math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
	end
	for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
		for i=k+1:N 
			<math>\displaystyle a_{ij}</math> -= <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 Rozwiązywanie układu równań metodą 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>;

<flash>file=eliminacjagaussa.swf</flash><div.thumbcaption>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 wykorzystujący pamięć podręczną?

Choć programista nie ma bezpośredniego wpływu na opisane powyżej działania systemu operacyjnego i hardware 'u (zarządzających, odpowiednio, pamięcią wirtualną i cache ), to przez właściwe projektowanie algorytmów -- a zwłaszcza: ich właściwą implementację -- może spowodować, że jego programy nie będą zmuszały komputera do nieracjonalnych zachowań. Jak się okazuje, całkiem łatwo nieświadomie tworzyć programy przeczące zasadom efektywnego wykorzystania hierarchii pamięci. Na potwierdzenie zacytujmy za Dongarrą \cite{Dongarra} klasyczny już przykład mnożenia dwóch macierzy.

W programie wykonujemy operację mnożenia dwóch macierzy $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

 [ijk]
/* 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:

 [ikj]
/* 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$:

 [bikj(16)]
/* 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. (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 $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$,

$\left(\begin{array}{ccc}{a}_{11}& \cdots & a_{1m}\\ ⋮& & ⋮\\ a_{n1}& \cdots & a_{nm}\end{array}\right).$

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:

• 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 $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. 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. \cite{clapack-howto}).

Przypomnijmy przy okazji, że ze względu na konstrukcję cache 'a spotykaną np. w procesorach Intela i AMD, czasem warto stosować tzw. array padding w przypadku, gdy mamy do czynienia z macierzami o wymiarach będących dużą potęgą dwójki, zob. Rozdział Uzupelnic: sec:cache:example .