MN06

Pamięć hierarchiczna komputerów a algorytmy numeryczne

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 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. 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). Oto jakie wyniki uzyskaliśmy dla obliczeń w arytmetyce podwójnej precyzji double na maszynie z procesorem AMD Duron i zegarem 1.1 GHz:

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 */
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.

Algorytm ikj

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

 
/* 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ę, zmniejszając liczbę cache misses!

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

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.

Reprezentacja macierzy gęstych

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

Dlatego 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) (pamiętajmy wszak, że 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 pokazujemy przykładowy fragment kodu realizującego opisaną wyżej ideę.

#define N 10
#define IJ(i,j,n) ((i)-1+((j)-1)*(n))
#include <stdio.h>

int main(void)
{
double *matrix, *ptr;
int i,j;

matrix = (double *)malloc(N*N*sizeof(double));

ptr = matrix;

/* staramy się jak naczęściej odwoływać się do elementów macierzy kolumnowo */

for (j=1; j<=N; j++)
	for (i=1; i<=N; i++)
	{
*ptr = i+2*j; 	/* przypisanie wartości elementowi macierzy */ 
		ptr++;		/* przejście do kolejnego elementu */
	}
		
/* jeszcze raz to samo, ale w inny, bardziej czytelny, choć mniej optymalny
sposób */

for (j=1; j<=N; j++)
	for (i=1; i<=N; i++)
	{
		matrix[IJ(i,j,N)] = i+2*j;
	}
		
/* dla wydruku, odwołujemy się do elementów macierzy wierszowo */
for (i=1; i<=N; i++)
{
	for (j=1; j<=N; j++)
		fprintf(stderr,"%5.2g ", matrix[IJ(i,j,N)]);
	fprintf(stderr,"\n");
}

free(matrix);		
return(0);
}

Dodajmy, że opisane podejście nie jest niczym nowym w praktyce numerycznej i jest stosowane w wielu dobrych bibliotekach numerycznych.

==Włączenie zewnętrznej biblioteki fortranowskiej do programu==

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.

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:

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

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:

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

No to wykorzystajmy naszą funkcję:

 
/* Wykorzystanie funkcji DNRM2 fortranowskiej biblioteki BLAS w programie w C*/
#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: %e\n", dnrm2_(&n, x, &incx));

	return(0);
}

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

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

 
/*Wykorzystanie funkcji fortranowskiej operującej na macierzy. 
Zwróćmy uwagę na sposób użycia argumentu macierzowego*/
#include <stdio.h>
#include <stdlib.h>
int fillmatrix_(int *, int *, double *); 

int main()
{
	int MM, NN, i, j;
	double *A;
	
   	MM = 2; NN = 5;
	
	A = (double *)malloc(MM*NN*sizeof(double));

   	fillmatrix_( &MM, &NN, A );

	printf("\nKolejne elementy wektora A:\n\n");   
   	for ( i = 0; i < NN*MM ; i++ ){
		printf("%e\n", A[i] );
	}

	printf("\nWektor A zinterpretowany jako macierz:\n\n");   
	for ( j = 0 ; j < MM ; j++ )
	{
   		for ( i = 0; i < NN ; i++ )
   			printf("%e ", A[i*MM+j] );
		printf("\n");
   	}

   	free( A );

   	return(0);
}

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 oraz LAPACK Linear Algebra PACKage . 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, . 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 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;

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

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

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

Programy korzystające z BLASów i LAPACKa 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