MN04

BLAS, LAPACK i ATLAS=

W zadaniach dotyczących macierzy gęstych, warto skorzystać z klasycznego tandemu bibliotek: BLAS (Basic Linear Algebra Subprograms ) \cite{BLAS-home-page} oraz LAPACK (Linear Algebra PACKage ) \cite{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, \cite{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 \cite{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 \cite{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

Mnożenie macierz-wektor w BLAS

Zacznijmy od prostej procedury BLAS Level 2, jaką jest mnożenie macierzy przez wektor. Wykorzystamy tzw. wzorcową implementację BLASów (niezoptymalizowaną) dostarczaną z dystrybucją np. Red Hat Linux. Jest to biblioteka funkcji fortranowskich.

Naszym zadaniem jest wykonanie operacji

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

 
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 dające wygodny dostęp do macierzy niezależny od jej wewnętrznej reprezentacji, a także 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:

#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 = 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. 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.