2007.10_Biblioteka funkcji numerycznych LAPACK_[Programowanie].pdf

dla programistów

Obliczenia numeryczne

numerycznych LAPACK

Marek Sawerwain

Obliczenia numeryczne to dziedzina, w której komputery odgrywają najważniejszą rolę, ponieważ bez

maszyn liczących, nie można przeprowadzać tysięcy, milionów, a nawet miliardów obliczeń na sekundę.

Jednak, co warto wiedzieć, podczas jakichkolwiek operacji na liczbach z przecinkiem, maszyny

popełniają błędy.

mać błędne wyniki, dalekie od prawdziwych lub

w najlepszym przypadku znacząco zniekształcone.

Ten problem wynika z ograniczonej precyzji, liczb,

na których wykonujemy obliczenia. Stosowanie rozwią-

zań, gdzie precyzja zapisu jest bardzo duża, np: tysiąc

miejsc po przecinku jest możliwe. Jednakże są to zawsze

rozwiązania programowe o niskiej wydajności w po-

równaniu do obliczeń, jakie procesor wykonuje na typach

danych zaimplementowanych sprzętowo.

Dlatego też z czasem opracowano algorytmy, któ-

re pomimo ograniczonej precyzji są dokładne oraz szyb-

kie. Jednym z pakietów, udostępnionych razem z ko-

dem źródłowym jest LAPACK, czyli zbiór procedur nu-

merycznych, napisanych w języku Fortran 77 . Język oraz

jego współczesne odmiany, jak np: Fortran 95 , to języki

naturalnie uniwersalne, ale stosowane głównie do obli-

czeń numerycznych.

Jednakże, dzięki pakietowi GCC , nie ma przeciwwska-

zań do wykorzystywania dostępnych funkcji, nie tylko

w Fortranie, ale również w języku C czy C++. Jest to uza-

sadnione tym, że darmowe biblioteki GUI, jak choćby QT

czy GTK+ nie mają portu języka Fortran . Oznacza to, że

aplikację użytkownika trzeba napisać z C/C++, ale pro-

cedury numeryczne mogą być bez żadnych przeszkód

tworzone w Fortranie . W tym artykule chciałbym pokazać

sposób korzystania z pakietu LAPACK, właśnie z pozio-

mu języków C/C++.

Kompilacja pakietów

Pakiet LAPACK oraz pozostałe elementy, jak BLAS czy

ATLAS, mogą być dostępne w repozytoriach dystrybucji,

której używamy, a przykładem są dystrybucje Debian oraz

Ubuntu. Jednakże może się okazać, iż nie mamy łatwego

dostępu do wersji binarnych. W takim przypadku trzeba

przeprowadzić kompilację.

Wszystkie podstawowe pakiety, takie jak LAPACK,

BLAS oraz CBLAS , niestety nie mają skryptów conigure ,

a tylko zawierają prosty system do kompilacji, oparty o pli-

ki makeile . Jednak, wbrew pozorom kompilacja jest dość

prosta, a jedynie trudniej przeprowadzić kompilację pa-

kietu ATLAS , bowiem choć oferuje skrypt conigure, to nie

jest tworzony z pomocą pakietów Automake i Autoconf , ale

swoje zadanie spełnia bardzo dobrze. Dodatkowo trzeba

październik 2007

Biblioteka funkcji

J eśli źle zapiszemy nasze obliczenia, to możemy otrzy-

dla programistów

Obliczenia numeryczne

Listing 1. Oryginalny nagłówek procedury numerycznej cheev

na nową:

FC = gfortran

SUBROUTINE CHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,

RWORK,INFO )

*– LAPACK driver routine (version 3.1) –

Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.

November 2006

Scalar Arguments

CHARACTER JOBZ, UPLO

INTEGER INFO, LDA, LWORK, N

Array Arguments .

REAL RWORK( * ), W( * )

COMPLEX A( LDA, * ), WORK( * )

Pakiet CBLAS nie zastępuje pakietu BLAS ,

więc przed rozpoczęciem kompilacji trzeba

pamiętać, aby biblioteka libblas.a znajdowała

się już w systemie. Analogicznie wykonuje-

my, gdy wykorzystujemy ten pakiet w swo-

ich programach. Podczas konsolidacji uży-

wamy pakietów BLAS jak i CBLAS .

Przygotowania do kompilacji pakie-

tu LAPACK, są dość podobne jak w przy-

padku pakietu CBLAS . Z podkatalogu IN-

STALL , trzeba skopiować plik make.inc.LI-

NUX do katalogu głównego LAPACK'a

i zmienić jego nazwę na make.inc . Podob-

nie jak poprzednio, należy zmienić nazwę

kompilatora z g77 na gfortran w pliku ma-

ke.inc . Podczas kompilacji warto wydawać

polecenie make z parametrami, jeśli chce-

my uzyskać pakiet BLAS to wydajemy po-

lecenie:

przeprowadzać specjalne testy, a w zależno-

ści od wersji pakietu kompilacja przebiega w

różny sposób.

Najszybciej można wykonać kompilację

pakietu BLAS , przy czym pakiet LAPACK

ma własną kopię BLAS'a , więc nie jest to ko-

nieczne. Jeśli jednak chcemy dokonać kompi-

lacji, to w pierwszej kolejności musimy ścią-

gnąć archiwum pakietu BLAS . Będzie to plik

o nazwie blas.tgz (plik jest niewielki około

100kb), a jego rozpakowanie wykonamy za

pomocą polecenia:

oraz:

LOADER = g77

Zmienić trzeba nazwę kompilatora z g77 na

gfortran :

FORTRAN = gfortran

make blaslib

oraz:

kompilację pakietu LAPACK rozpoczniemy

od polecenia:

LOADER = gfortran

make lapacklib

tar zxvpf blas.tgz

I to wszystkie zmiany. Po wydaniu polecenia

make i po kilku chwilach zostanie utworzona

biblioteka pakietu BLAS .

W dość prosty sposób dokonujemy kom-

pilacji pakietu CBLAS , czyli specjalnego API

do korzystania z funkcji BLASA w języku C.

Warto z tej biblioteki korzystać, ponieważ uła-

twia używanie funkcji z pakietu BLAS w ję-

zyku C. Po dekompresji archiwum, nim wy-

damy polecenie make, należy wykonać jed-

ną czynność, tzn. wskazać plik z konigura-

cją, czyli nazwą kompilatora. Zakładamy, że

kompilujemy bibliotekę pod systemem Linux,

więc wykorzystujemy plik o nazwie Makei-

le.LINUX , więc musimy zmienić nazwę tego

pliku na Makeile .in lub po prostu wykonać je-

go kopię np:

System kompilacji składa się z dwóch pli-

ków Makeile oraz make.inc . Jeśli używamy

kompilatora GCC w wersji 3.X, możemy od

razu wydać polecenie make i po kilku chwi-

lach otrzymamy plik z gotową biblioteką

o nazwie blas_LINUX.a. Jednakże we wszyst-

kich nowych dystrybucjach systemu Linux

obecny jest kompilator w wersji 4.X. Nato-

miast jedną z wielu zmian jest zmiana nazwy

kompilatora Fortrana z g77 na gfortran. Z te-

go powodu musimy za pomocą jakiegokol-

wiek edytora tekstowego zmienić dwie linie

w pliku make.inc .

Obydwie linie przedstawiają się nastę-

pująco:

Nie warto wydawać polecenia make bez

parametrów, ponieważ po kompilacji roz-

poczną się testy, które trwają dość długo

i w większości przypadków nie są potrzebne

do przeprowadzania. Sam proces kompila-

cji zabiera sporo czasu choćby dlatego, że

w katalogu SRC znajdują się 1348 pliki źró-

dłowe w Fortranie .

Korzystanie z funkcji Fortrana

w języku C

Stosując kompilator GCC, możemy w dość

łatwy sposób korzystać z funkcji LAPACKA,

w projektach tworzonych w języku C. Istnie-

je nawet pakiet CLAPACK , który podobnie

Listing 2. Dodatkowo wartość parametru INFO

FORTRAN = g77

cp Makeile.LINUX Makeile.in

Listing 6. Nagłówek funkcji numerycznej snrm2

albo utworzyć link symboliczny:

int cheev ( char JOBZ , char UPLO ,

int N , complex * A , int LDA , loat

* W , complex * WORK , int LWORK , loat

* RWORK ) {

int INFO ;

cheev_ (& JOBZ , & UPLO , & N , A , & LDA ,

W , WORK , & LWORK , RWORK , & INFO );

return INFO ;

}

REAL FUNCTION SNRM2(N,X,INCX)

* . Scalar Arguments

INTEGER INCX,N

* .

* . Array Arguments

REAL X(*)

ln -s Makeile.LINUX Makeile.in

analogicznie jak w przypadku pakietu BLAS ,

trzeba jeszcze zmienić nazwę kompilatora

Fortrana, z nazwy:

FC = g77

www.lpmagazine.org

dla programistów

Obliczenia numeryczne

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

cierz to zwykły wskaźnik na ciąg wartości

complex, ale trzeba zwracać baczną uwa-

gę na organizację macierzy, a mianowi-

cie, czy macierz jest zorientowana wierszo-

wa czy kolumnowa. Bardzo często, w funk-

cjach i procedurach obowiązuje orientacja

kolumnowa.

Możemy utworzyć pomocniczą funk-

cję o nazwie cheev , zbudujemy ją w ten spo-

sób, iż część argumentów będziemy prze-

kazywać przez wartość, a tylko macierze

będą przekazywane za pomocą wskaźni-

ków. Dodatkowo wartość parametru INFO

będzie stanowić wynik działania funkcji

cheev . Listing 2

Wartość INFO jest bardzo ważna, ponie-

waż za pomocą tego parametru, LAPACK in-

formuje, czy obliczenia zostały wykonane po-

prawnie.

W pakietach LAPACK oraz BLAS , oprócz

podprogramów, obecne są także funkcje,

a przykładem może być funkcja snrm2 , ob-

liczająca euklidesową normę z podane-

go wektora. Mamy trzy argumenty, dwie

liczby całkowite oraz wektor składający się

z typu REAL.

Nagłówek tej funkcji w języku C przed-

stawia się następująco:

int N = 4 , NRHS = 1 ;

loat A []= { 1 , 2 , 3 , 3 ,

2 , 3 , 4 , 5 ,

2 , 2 , 3 , 4 ,

4 , 3 , 2 , 2 } ;

int LDA = 4 ; int IPIV [ 4 ]= { 0 , 0 , 0 , 0 } ;

loat B [ 4 ]= { 2 , 3 , 3 , 2 } ;

int LDB = 4 , INFO ;

void print_vec_loat ( int x , loat * tab ) {

int i , j ;

for ( i = 0 ; i < x ; i ++) {

printf ( "%2.2f " , *( tab + i ));

}

printf ( " \n " );

}

transpose_loat ( 4 , 4 , & A [ 0 ]);

INFO = sgesv ( N , NRHS , & A [ 0 ] , LDA , & IPIV [ 0 ] , & B [ 0 ] , LDA );

if ( INFO == 0 ) print_vec_loat ( 4 , & B [ 0 ]);

jak CBLAS oferuje dość wygodne API. Jed-

nak ten pakiet dotyczy wersji 3.0, a najnow-

sza to 3.1.1, ponadto jak się za chwilę okaże,

uzyskanie dostępu do poszczególnych funk-

cji jest bardzo proste i można zrealizować sa-

modzielnie.

Załóżmy, że chcemy skorzystać z funkcji

o nazwie cheev . Zadaniem jej jest wyznacze-

nie wartości własnych oraz wektorów wła-

snych zespolonej macierzy hermitowskiej. Li-

sting 1 pokazuje nagłówek funkcji cheev , jaki

można znaleźć w pliku źródłowym cheev.f.

Choć nazwa funkcji, czy raczej procedu-

ry, jest zapisana za pomocą dużych liter, to

w rzeczywistości zostanie zakodowana ja-

ko cheev (ze znakiem podkreślenia na koń-

cu). Typy parametrów, jakie są przekazy-

wane, mają, co jest dla nas bardzo korzyst-

ne, bezpośrednie odniesienie do języka C,

np: typ REAL w Fortranie to typ loat . Zbior-

cze podsumowanie, wszystkich podstawo-

wych typów zawiera tabela pt: Podstawowe

typy danych w Fortranie i ich odpowied-

niki w C.

Charakterystycznym elementem jest to,

iż wszystkie parametry są przekazywane ja-

ko pojedyncze wskaźniki. Nie jest ważne,

czy jest to pojedyncza zmienna, czy też ma-

cierz. Te informacje wystarczą, aby popraw-

nie zapisać prototyp funkcji cheev o następu-

jącej postaci:

W ten sposób uzyskaliśmy gotową do wyko-

rzystania funkcję w języku C. Jak widać uzy-

skanie dostępu sprowadza się tylko do od-

powiedniego oprogramowania argumen-

tów. Wygodne jest także przekazywanie ma-

cierzy, np: parametr A w funkcji cheev . Ma-

Listing 7. Wyznaczenie pierwiastka z macierzy hermitowskiej

complex * A , * WORK ;

loat * W , * RWORK , entropy ;

int N = 2 , INFO , LDA , LWORK ;

LDA = N ; LWORK = 2 * N - 1 ;

W = malloc ( N * sizeof ( loat ));

WORK = malloc ( LWORK * sizeof ( complex ));

RWORK = malloc (( 3 * N - 2 )* sizeof ( loat ));

A = malloc ( N * N * sizeof ( complex ));

INFO = cheev ( 'V', 'U', N , A , LDA , W , WORK , LWORK , RWORK );

_Complex sqrt_cc , in_data ;

diag = create_matrix ( N , N );

for ( i = 0 ; i < N ; i ++) {

if (*( W + i ) > 0 ) set_cell_at_matrix_direct ( diag , i , i , sqrt (*( W + i )) , 0 );

if (*( W + i ) == 0 ) set_cell_at_matrix_direct ( diag , i , i , 0 , 0 );

if (*( W + i ) < 0 ) {

__real__ in_data = *( W + i );

__imag__ in_data = 0 ;

sqrt_cc = csqrtf ( in_data );

set_cell_at_matrix_direct ( diag , i , i , crealf ( sqrt_cc ) , cimagf ( sqrt_cc ));

}

transpose_matrix ( a_mat );

tmp_mat = create_matrix ( a_mat -> rows , a_mat -> cols );

tmp1 = create_matrix ( a_mat -> rows , a_mat -> cols );

mul_matrix ( a_mat , diag , tmp1 );

transpose_matrix ( a_mat );

mul_matrix ( tmp1 , a_mat , tmp_mat );

extern void cheev_ (char *JOBZ,

char *UPLO, int *N, complex *A,

int *LDA, loat *W, complex

*WORK, int *LWORK, loat *RWORK,

int *INFO);

październik 2007

dla programistów

Obliczenia numeryczne

extern loat snrm2_(int *N, loat *X,

int *INCX);

znajduje się tam zero, oznacza to, że LAPACK

wykonał wszystkie obliczenia poprawnie.

Jeśli mamy wartość większą od zera, to na

pozycji (i, i) rozkładu LU macierzy A (A =

P*L*U) pojawiło się zero, co niestety oznacza,

iż mamy macierz osobliwą.

Również wartość ujemna oznacza błąd,

bowiem na pozycji i jeśli pozbędziemy się zna-

ku minus, znajduje się niepoprawna wartość.

Zakładając, że wszystkie obliczenia się powio-

dły, rozwiązanie naszego układu znajduje się

w zmiennej B, natomiast wyświetleniem wyni-

ku zajmuje się funkcja print_vec_loat.

Pakiet ATLAS

natomiast wersja ugrzeczniona , bez wskaźni-

ków w typach skalarnych, przedstawia się

następująco:

Pakiet ATLAS, czyli Automatically Tuned Li-

near Algebra Software, to wyspecjalizowany

pakiet, którego zadaniem jest dostarczenie

bardzo wydajnych procedur do operacji na

macierzach i wektorach. Pakiet LAPACK do

realizacji podstawowych operacji na macier-

zach, stosuje funkcje pakietu BLAS. Pod-

miana najważniejszych funkcji tego pakie-

tu na wydajniejsze, powoduje zwiększa-

nie wydajności pozostałych funkcji LAPAC-

K'a, a przykładem może być funkcja matmul

(mnożenie macierzy) zastąpienie wydajniej-

szą wersją, powoduje wzrost wydajności

wszystkich funkcji, należących do tzw. po-

ziomu trzeciego w pakiecie BLAS.

Zysk na wydajności jest bardzo duży,

a osiągnięty został poprzez dopasowanie

poszczególnych funkcji do specyicznych wła-

ściwości procesorów. Wykorzystywane są

instrukcje SSE, które pozwalają nawet na

czterokrotne przyspieszenie wykonywania

niektórych operacji na macierzach. Wyko-

rzystywane są także instrukcje 3DNow, ofe-

rowane przez procesory AMD. Natomiast

ATLAS oferuje wsparcie innych procesorów,

np.: wykorzystywane są instrukcje Altivec

oraz specyiczne własności procesorów Po-

werPC czy UltraSparc. Warto wykonać sa-

modzielną kompilację tego pakietu w ra-

mach swojego systemu, ponieważ jeśli ma-

my procesor z większą ilości pamięci cache,

to ATLAS dopasuje się do możliwości na-

szego procesora i w efekcie wydajność mo-

że być jeszcze większa.

loat snrm2(int N, loat *X, int INCX) {

return snrm2_(&N, X, &INCX);

}

Rozwiązywanie układu równań

liniowych

Rozwiązywanie układu równań liniowych

wykonuje funkcja sgesv . W zależności od spe-

cyicznej postaci macierzy, LAPACK oferu-

je inne, bardziej efektywne metody rozwią-

zywania układów równań typu AX=B. Jed-

nak w przypadku, gdy mamy pełną macierz,

to stosujemy funkcję sgesv , która rozwiązu-

je najbardziej ogólny przypadek. Nasz układ

równań liniowych jest następujący:

Obliczanie pierwiastka z macierzy

Rozwiązanie układu równań liniowych spro-

wadziło się do wywołania jednej funkcji. LA-

PACK, która została utworzona po to, aby

typowe zadania numeryczne rozwiązywać

za pomocą pojedynczych funkcji. Jednak-

że nie jest tak zawsze, np: jeśli chcemy obli-

czyć pierwiastek z macierzy (jeśli A jest ma-

cierzą, to pierwiastek określamy jako q=Sqr-

t(A), podobnie jak dla liczb q*q=A), to nie

znajdziemy jednej funkcji, która ułatwi roz-

wiązanie.

Załóżmy, że chcemy obliczyć pierwia-

stek z hermitowskiej macierzy o wartościach

zespolonych. W pierwszej kolejności musi-

my wyznaczyć wartości własne tej macierzy

oraz macierz wektorów własnych. Bowiem

taka dekompozycja pozwala zapisać macierz

A jako (V – macierz zbudowana z wekto-

rów własnych, d – macierz zbudowana z war

tości własnych umieszczonych na głównej

przekątnej, oraz V – to transponowana ma-

cierz V):

• 1a + 2b + 3c + 3d = 2

• 2a + 3b + 4c + 5d = 3

• 2a + 2b + 3c + 4d = 3

• 4a + 3b + 2c + 2d = 2

Nim wywołamy właściwą funkcję, to trzeba

zadeklarować odpowiednie zmienne, opisu-

jące nasz układ. Przez N oznaczamy wymiar

macierzy. W naszym przypadku wartość

N to 4. Kolejna zmienna NRHS opisuje ile ko-

lumn znajduje się w macierzy B, czyli wy-

razów wolnych.; w naszym przypadku ma-

my tylko jedną kolumnę. Tworzymy także

zmienną LDA , gdzie znajdzie się liczba wier-

szy macierzy A. Podobną rolę pełni zmien-

na LDB , ale jest to liczba wierszy macierzy B

(wyrazy wolne). Deklarujemy także jeszcze

tablicę IPIV , w której znajdą się informacje

o poprzestawianych wierszach.

Po deklaracji wszystkich zmiennych moż-

na już wywołać sgesv , która rozwiąże to rów-

nanie. Niestety, okaże się, iż otrzymaliśmy wy-

niki, które nie pasują do poprawnego rozwią-

zania. Powód jest bardzo prozaiczny – sgesv

wymaga, aby przekazać macierz w orienta-

cji kolumnowej. Możemy przepisać deklarację

zmiennych lub wykonać transpozycję macie-

rzy A i właśnie tym zajmuje się funkcja trans-

pose_loat, której wywołanie następuje przed

funkcją sgesv. Listing 6 przedstawia najważ-

niejsze fragmenty kodu rozwiązującego nasz

układ równań (pełny kod źródłowy przykła-

dów znajduje się na płycie DVD).

Po wykonaniu funkcji sgesv , warto spraw-

dzić wartość wynikową tej funkcji, czyli wspo-

minaną już wcześniej wartość INFO. Gdy

A = V*d*V'

wolną potęgę macierzy, a jak łatwo zgadnąć,

będzie to bez porównania szybszy sposób niż

kolejne mnożenia macierzy.

obliczenie pierwiastka polega na modyika-

cji wartości własnych, a dokładniej na poli-

czeniu pierwiastka z każdej wartości własnej,

umieszczonej na głównej przekątnej:

A = V * (d^n) * V'

A = V * SQRT(d) * V'

Do rozłożenia macierzy hermitowskiej po-

trzebna jest funkcja cheev , która działa na licz-

bach zespolonych complex . Zakładamy, że bę-

dziemy stosować pojedynczą precyzję.

w ten sposób możemy obliczyć również do-

Tabela 1. Podstawowe typy danych w Fortranie i ich odpowiedniki w C

Typ Fortrana

Typ języka C

CHARACTER

char

INTEGER

int

COMPLEX

struct { loat r, i; } complex;

DOUBLECOMPLEX

struct { double r, i; } doublecomplex;

REAL

loat

DOUBLEREAL

double

LOGICAL

long int

www.lpmagazine.org

dla programistów

Obliczenia numeryczne

O autorze

Autor zajmuje się tworzeniem oprogramo-

wania dla WIN32 i Linuksa, jego zaintere-

sowania to teoria języków programowania

oraz dobra literatura.

Kontakt z autorem: ay@linux.com.pl.

__real__ in_data = *(W+i);

__imag__ in_data = 0;

sqrt_cc = csqrtf(in_data).

Rysunek 1. Podstawowa strona, dotycząca pakietu LAPACK

Podobnie jak w poprzednim przykła-

dzie, trzeba zadbać o parametry. Zmienne ty-

pu N, LDA mają takie same znaczenie jak po-

przednio, a główna macierz przechowywana

jest w zmiennej A, wartości własne po oblicze-

niu będą przechowywane w zmiennej W, na-

tomiast macierz zbudowana z wektorów wła-

snych znajdzie się na miejscu oryginalnej ma-

cierzy A.

Do poprawnej pracy potrzebna jest dodat-

kowa tablica WORK oraz RWORK . Jak widać

na Listingu 4 część zmiennych jest deklarowa-

na jako wskaźniki, a także przydzielana jest pa-

mięć, za pomocą polecenia malloc. Nie prze-

szkadza to funkcji cheev w poprawnym dzia-

łaniu, ponieważ samo wywołanie funkcji jest

następujące:

Z obliczonych wartości budujemy macierz,

gdzie na przekątnej znajdą się potrzebne

wartości. Ostatnia część Listingu 4 to wyko-

nanie ostatecznych mnożeń. Obliczenia są re-

alizowane od końca, pierwsze mnożenie to

d*V', otrzymany wynik jest mnożony przez

V. Dlatego transpozycja jest wykonywana

w pierwszej kolejności, a po wykonaniu

mnożenia wykonujemy transpozycję raz

jeszcze, aby dokonać drugiego mnożenia.

Historia pakietu LAPACK

Pakiet LAPACK ma dość długą historię; mo-

żna nawet powiedzieć, iż niespotykanie dłu-

gą, biorąc pod uwagę szybki rozwój informa-

tyki. Pierwsza wersja pojawiła się w 1992 r.,

natomiast wersja 2.0 w 1994 r., potem na-

stąpiła przerwa aż do 1999 r., kiedy wyda-

no wersję 3.0. Po tym czasie znów nastąpi-

ła przerwa, bowiem na następną wersję 3.1

trzeba było czekać aż do 2006 r. Jednak na

początku 2007 r. pojawiło się kolejne uaktu-

alnienie i obecna wersja pakietu 3.1.1.

LAPACK to około 800 tys. linii kodu

w Fortranie 77, gdzie funkcje zostały podzie-

lone na cztery podstawowe typy danych:

real, complex, double precision, double com-

plex. Oznacza to, iż wiele funkcji jest imple-

mentowanych czterokrotnie, a każda imple-

mentacja bierze pod uwagę specyiczne

własności typu na jakim pracuje.

Wraz z funkcjami numerycznymi, w pa-

kiecie LAPACK znajdziemy kody źródłowe

pakietu BLAS. Pakiet, o pełnej nazwie Ba-

sic Linear Algebra Subprograms, zawiera

pomocnicze funkcje, które wykonują podsta-

wowe operacje na macierzach, a wydajność

funkcji LAPACK, ściśle zależy od ich wydaj-

ności. Zalecane jest, aby używać specjal-

nie dedykowanych pakietów, np.: do kon-

kretnych procesorów AMD czy INTEL. Pa-

kiet BLAS ma również imponującą historię,

warto przejrzeć kod źródłowy funkcji BLA-

S'a, aby przekonać się, iż kod powstał w

drugiej połowie 70. lat , a pomimo upływu 30

lat nadal oferuje bardzo wygodne i użytecz-

ne podprogramy i funkcje.

Podsumowanie

Ilość funkcji oferowana przez pakiet LAPACK

jest bardzo duża, jednak opisane w tym artyku-

le przykłady pokazują, w jaki sposób możemy

korzystać z tej znakomitej biblioteki, aby wy-

konywać we własnych programach obliczenia

numeryczne. Mimo, że powstają nowe biblio-

teki do obliczeń, to w zakresie algebry liniowej,

czyli uogólniając, w zastosowaniach związa-

nymi z obliczeniami na macierzach, LAPACK

nadal pozostaje najlepszą biblioteką. Warto,

wypróbować ten pakiet we własnych projek-

tach, a pewnością w ten sposób zaoszczędzi-

my sporo czasu.

INFO = cheev('V', 'U', N, A, LDA, W,

WORK, LWORK, RWORK);

Pierwszy argument 'V' oznacza, że chcemy

obliczyć wektory oraz wartości własne. Dru-

gi U oznacza, iż macierz jest tak zapisana, że

w pierwszej kolejności podana została część

macierzy, leżąca nad główną przekątną. War-

tość ta wskazuje, czy po obliczeniu tylko warto-

ści własnych ma zostać zniszczona górna bądź

dolna część macierzy. Takie podejście oznacza

mniejsze zapotrzebowanie na pamięć, choć jak

widać, część macierzy jest tracona, albo cała

macierz jest zastępowana przez wektory własne.

Po wykonaniu obliczeń możemy przy-

stąpić do dalszych działań, a są nimi oblicze-

nie pierwiastków z wartości własnych. Ponie-

waż mogą pojawić się ujemne wartości, musi

być przygotowana na obliczanie pierwiastków

z ujemnych liczb, dlatego używamy kompila-

tora GCC, i wykorzystamy wbudowany typ

_Complex oraz specjalne operatory __real__

oraz __imag__ do odczytu części rzeczywi-

stej i urojonej. Obliczanie pierwiastka z liczby

ujemnej następuje przy użyciu funkcji csqrtf,

która stanowi rozszerzenie dostępne w kom-

pilatorze GCC:

W Sieci

• Podstawowa strona, dotycząca

pakietu LAPACK:

http://www.netlib.org/lapack/

• Port pakietu LAPACK do języka C:

http://www.netlib.org/clapack/

• Podstawowa strona, dotycząca

pakietu BLAS i CBLAS:

http://www.netlib.org/blas/

• Podstawowa strona, dotycząca

pakietu ATLAS:

http://math-atlas.sourceforge.net/

• Przykłady w Fortranie, rozwiązujące

różne zagadnienia numeryczne,

przy równoczesnym wykorzystaniu

funkcji pakietu LAPACK:

http://www.nag.com/lapack-ex/

październik 2007

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika: