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:

1  

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

1. 1

(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:

1  

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

1. 1

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:

1  [Algorytm eliminacji Gaussa]

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

1. 1

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ą).

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

1. 1

Ł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ć {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.

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

1. 1

Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest już rozkład LU wykonany z wyborem w kolumnie.

1  [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$;

1. 1

{}{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

1  [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];

1. 1

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:

1  [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];

1. 1

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

1  [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];

1. 1

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

===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 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 f2c lub jej modyfikacji na użytek kompilatora GCC, biblioteki gfortran .

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 liblapack.a , 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 DNRM2 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 := sqrt( x'*x )

i tak dalej...

Tak więc nasza funkcja obliczająca normę wektora ma trzy argumenty: N -- długość wektora ( INTEGER ), X -- wektor, którego długość chcemy obliczyć (tablica liczb DOUBLE PRECISION ) oraz tajemniczy dodatkowy parametr INCX typu INTEGER -- jest to wartość skoku, określająca co który element wektora uwzględnić w obliczaniu normy: aby policzyć normę całego wektora, bierzemy INCX równe 1. Używając zapisu Octave, DNRM2 oblicza po prostu

norm( X(1:INCX:N) )

Kod obiektowy tej funkcji znajduje się już w bibliotece BLAS, zawartej w pliku libblas.a . 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 gcc / gfortran , nazwą funkcji do wykorzystania w C będzie dnrm2_ (tak! małymi literami i z przyrostkiem " _ ").

Jeśli zaś chodzi o argumenty, to zapewne co do drugiego nie będziemy mieli wątpliwości: jako wektor X przekażemy -- naturalnie -- wskaźnik do tablicy X (typu double ), czyli po prostu: jej nazwę. Co z pozostałymi argumentami? Okazuje się, że reguła jest niezmienna:

Każdy argument funkcji fortranowskiej zastępujemy wskaźnikiem do odpowiedniego typu:

Uzupelnij tytul

Fortran 77 || C
INTEGER || int
REAL	float
DOUBLE PRECISION	double
COMPLEX	struct { float Re, Im; }
DOUBLE COMPLEX	struct { double Re, Im; }
CHARACTER	char

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

A więc pierwszym i trzecim argumentem funkcji dnrm2_ będą wskaźniki do int . Ponieważ funkcja DNRM2 zwraca w wyniku liczbę podwójnej precyzji, to ostatecznie prototyp naszej funkcji w C będzie następujący:

1  

double dnrm2_(int* N, double* X, int* INCX);

1. 1

No to wykorzystajmy naszą funkcję:

1  [Wykorzystanie funkcji DNRM2 fortranowskiej biblioteki BLAS w programie w C]

1. include <stdio.h>

double dnrm2_(int*,double*,int*);

int main(void) { int n, incx=1; double x[3]= {0,1,2};

n = 3; printf("Norma podanego wektora: return(0); }

1. 1

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ę: gfortran , 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 $M\times N$:

[Funkcja w Fortranie 77, wypełniająca macierz $M\times N$]

SUBROUTINE FILLMATRIX ( M, N, MATRIX ) INTEGER M,N DOUBLE PRECISION MATRIX(M, N)

DO 10 I=1,M DO 20 J=1,N MATRIX(I,J) = 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 $M=2$, $N=5$, będzie macierz

Parser nie mógł rozpoznać (nieznana funkcja „\lstF”): {\displaystyle \displaystyle \lstF{MATRIX} = \begin{pmatrix} 3 & 5 & 7 & 9 & 11 \\ 4 & 6 & 8 & 10 & 12 \end{pmatrix} }

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 fillmatrix_ umieszczamy w osobnym pliku nagłówkowym ffortran.h , gdzie mamy zamiar kolekcjonować prototypy w C lokalnie wykorzystywanych funkcji fortranowskich):

1  [Wykorzystanie funkcji fortranowskiej

operującej na macierzy. Zwróćmy uwagę na sposób użycia argumentu macierzowego]

1. include <stdio.h>
2. include <stdlib.h>

int fillmatrix_(int *, int *, double *);

int main() { int MM, NN, i, j; double *A;

  	MM = 2; NN = 5;

A = (double *)malloc(MM*NN*sizeof(double));

  	fillmatrix_( &MM, &NN, A );

printf(" elementy wektora A:");

  	for ( i = 0; i < NN*MM ; i++ ){

printf("}

printf(" A zinterpretowany jako macierz:"); for ( j = 0 ; j < MM ; j++ ) {

  		for ( i = 0; i < NN ; i++ )
  			printf("printf("");
  	}

  	free( A );

  	return(0);

}

1. 1

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

{{{3}}}

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.