MN04

Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki

Rozwiązywanie układów równań liniowych

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

---

uuu

uuu

Algorytm ikj

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

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

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

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

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

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:

No to wykorzystajmy naszą funkcję:

Zwróćmy uwagę na sposób kompilacji tego programu:

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

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

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

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:

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:

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