A GSL FÜGGVÉNYKÖNYVTÁR TESZTELÉSE

A GSL FÜGGVÉNYKÖNYVTÁR TESZTELÉSE

Tartalomjegyzék

Bevezetés  3

1.1. Mi az a GSL 3

1.2. GSL telepítése 3

1.3. GSL használata 4

Vektorok és mátrixok 6

2.1. Vektorok és mátrixok készítése 6

2.2. Műveletek 11

2.3. BLAS 12

2.4. Normák 14

Transzformációk  15

3.1. QR felbontás 15

3.2. Householder transzformáció 16

Sajátérték, sajátvektor 18

Interpolációk  19

1. Bevezetés

1.1. Mi az a GSL

A GSL a GNU Scientific Library rövidítése. Ez egy olyan gyűjtemény, mely numerikus számításokat tartalmaz programozók számára. A GSL sztenderd C nyelven íródott, ezért minden benne található szolgáltatást függvényeken keresztül érhetünk el.

A következő fontosabb matematikai elemek találhatók benne:

Komplex számok
Lineáris algebra
Véletlen számok
Differenciál egyenletek
Numerikus differenciálás
Gyökkeresés
Vektorok és mátrixok
Sajátértékek
Interpoláció
Chebyshev egyenlőtlenség
stb

A GSL teljesen szabadon felhasználható és fejleszthető, így ha valaki kedvet kap a gyűjtemény továbbfejlesztéséhez az nyugodtan megteheti azt.

A bemutatást, tesztelést a GSL 1.6-os verziójával végezzük el.

A GSL könyvtár számos beépített funkcióval rendelkezik, úgyhogy mi nem azon leszünk, hogy egyes funkcióit felhasználva konstruáljunk olyan más funkciókat melyek szintén beépítve megtalálhatóak benne, hanem bemutatjuk azok használatát és teszteljük működésüket.

Persze nem teljes körű bemutatást végzünk csak belecsipegetünk a GSL nyújtotta szolgáltatásokba, aki komolyabban szeretne elmélyülni benne az számos leírást talál az interneten.

Jelen dolgozathoz az interneten fellelhető (www.gnu.org/software/gsl/) GSL-hez készült manuált használtuk. A matematikai hátérről pedig az olvasónak Virágh János: Numerikus matematika című könyvét ajánlhatjuk.

1.2. GSL telepítése

A telepítés LINUX alá könnyen véghezvihető. A GSL tartalmaz egy install nevű fájlt, melyben részletes le van írva a telepítés mibenléte. Az egész telepítés 4 konzolparanccsal elvégezhető:

./configure

make

make check

make install

Fontos azonban, hogy a megfelelő programok telepítve legyenek, főleg melyek a fordításhoz szükségesek: gcc, automake

1.3. GSL használata

A telepítés sikeres elvégzése után már semmi sem akadályozhat abban, hogy elkezdjünk komolyabban belemerülni a GSL nyújtotta szolgáltatásokba. Azonban előtte érdemes még szólni arról, hogy hogyan is használhatjuk fel a GSL-t a programjainknál.

Programok készítéséhez elegendő egy Nedit nevű programocska és a gcc C fordító. Mindkettőt a linux telepítésekor vagy utólag is feltehetjük gépünkre. Mi ezekkel fogunk dolgozni a GSL tesztelésekor.

Legelőször minden programunkban benne kell legyen a következő sorok valamelyike attól függően mit akarunk használni a programunkban:

#include <gsl/gsl_math.h>  // az alap típusokat tartalmazza

#include <gsl/gsl_vector.h>  // a vektorműveletekhez

#include <gsl/gsl_matrix.h>  // a mátrixműveletekhez

#include <gsl/gsl_blas.h>  // a BLAS függvéynek

#include <gsl/gsl_linalg.h>  // lineáris algebra

#include <gsl/gsl_eigen.h>  // sajátérték, sajátvektor

#include <gsl/gsl_interp.h>  // interpoláció

Ezzel máris használhatjuk a beépített függvényeket.

Egy másik fontos dolog a fordítás. A következő paranccsal fordíthatjuk le a programunkat:

gcc -static [forrásfájl].c -o [forrásfájl] -lgsl -lgslcblas -lm

A forrásfájl az általunk készített fájl neve kell hogy legyen. Sikeres fordítás esetén kapunk egy a forrásfájl nevével megegyező fájlt. Ez maga a lefordított program. A futtatáshoz csak be kell gépelnünk:

/[forrásfájl]

Ezzel készen is vagyunk arra, hogy elkezdjük tesztelni a GSL által nyújtott szolgáltatásokat.

2. Vektorok és mátrixok

2.1. Vektorok és mátrixok készítése

Vektorok és mátrixok készítéséhez a GSL tartalmaz két beépített adattípust. A vektorokhoz gsl_vector-t a mátrixokhoz gsl_matrix-t. Nézzünk egy példát vektor létrehozásához és adatokkal való feltöltéséhez valamin kiiratáshoz:

01. #include <stdio.h>

02. #include <gsl/gsl_vector.h>

04. void print_vector(gsl_vector* v)

14. int main(void)

A példából látható, hogy létrehozunk egy v vektort a gsl_vector_alloc(elemek_száma) segítségével, majd a gsl_vector_set(vektor,i,x) függvénnyel feltöltjük bizonyos értékekkel. Ha a 19. sorban alloc helyett calloc-ot írunk akkor a vektor 0 kezdőértékkel lesz feltöltve. A gsl_vector_free(vektor) felszabadítja a vektor által lefoglalt memóriát.

Speciális vektorok készítéséhez használhatók a következő függvények:

gsl_vector_set_zero(vektor)  – a vektor nullvektor lesz
gsl_vector_set_all(vektor,x) – a vektor minden eleme x lesz
gsl_vector_set_basis(vektor,i) – a vektor bázisvektor lesz, a i. elem lesz 1

Egy v vektor méretét a v->size attribútumával érhetjük el.

Most nézzük meg mi történik akkor, ha 09. sorban rosszul adjuk meg a vektor méretét. A következő üzenetet kapjuk:

gsl: vector_source.c:45: ERROR: index out of range
Default GSL error handler invoked.

A GSL minden hiba esetén azonnal megszakítja a program futtatását és megadja a hiba pontos helyét és a hiba rövid leírását. A GSL lehetőséget biztosít arra, hogy saját hibakezelő függvényeket hozzunk létre, melyek hiba esetén meghívódnak.

Mátrixokat is hasonlóan hozhatunk létre mint vektorokat. A létrehozáshoz használhatjuk a gsl_matrix_alloc(n,m) vagy a gsl_matrix_calloc(n.m) függvények valamelyikét. A mátrix megszüntetése a gsl_matrix_free(mátrix)-al lehetséges.

Mátrix elemeinek értéket a gsl_matrix_set(mátrix,i,j,x)-el adhatunk, lekérni pedig a gsl_matrix_get(mátrix,i,j) függvénnyel tudunk. Már itt megemlítjük, hogy hibás indexelés esetén a GSL megszakítja a program futását és egy hibajelzéssel leáll.

Különböző speciális mátrixok létrehozását a következő függvényekkel tehetjük meg:

gsl_matrix_set_all(mátrix,x) - a mátrix minden eleme x lesz
gsl_matrix_set_zero(mátrix) - a mátrix a nullmátrix lesz
gsl_matrix_set_identity(mátrix) - a mátrix az egységmátrix lesz

Egy m mátrix méretét az m->size1 és az m->size2 adattag adja meg.

Most nézzünk meg néhány különleges vektor illetve mátrixműveleteket a kódrészletekkel együtt.

Legyen v=[1 2 3 4 5 6 7 8 9], ez lesz a bemenetünk

KÓDRÉSZLET:
gsl_vector_swap_elements(v,0,1);
KIMENET:
v=[2 1 3 4 5 6 7 8 9]

A gls_vector_swap_elements(vektor,i,j) függvénnyel az i. elemet felcserélhetjük a j. elemmel.

KÓDRÉSZLET:
gsl_vector_reverse(v);
KIMENET:
v=[9 8 7 6 5 4 3 2 1]

A gsl_vector_reverse(vektor) függvény megfordítja a vektorunkat.

KÓDRÉSZLET:
gsl_vector_view view,view2;
gsl_vector* vv;
view = gsl_vector_subvector(v,1,4);
view2 = gsl_vector_subvector_with_stride(v,1,2,3);
vv = &(view.vector);
gsl_set_vector(vv,0,0);
KIMENET:
view=[2 3 4 5]
view2=[2 4 6]
vv=[2 3 4 5]
v=[0 2 3 4 5 6 7 8 9]

A gsl_vector_subvector(vektor,i,n) és a gsl_vector_subvector_with_stride(vektor,i,o,n) függvények mind az eredeti v vektorból részvektorokat képeznek. Azonban ezeket mind az úgynevezett gsl_vector_view típusba teszi meg. Ezt átkonvertálni gsl_vector típusba a példában látható módon lehetséges. Ezek után az átkonvertált részvektorral ugyanúgy végezhetőek műveletek, azonban arra vigyázni kell, hogy a módosítások az eredeti vektort módosítják. Ezt kiküszöbölhetjük, ha a részvektort átmásoljuk egy másik vektorba, de erről majd a következő részben lesz szó.

Legyen m=[[1 2 3] [4 5 6] [7 8 9]]

KÓDRÉSZLET:
gsl_vector* v1,v2;
v1 = gsl_vector_alloc(3);
v2 = gsl_vector_alloc(3);

gsl_matrix_get_row(v1,m,0);
gsl_matrix_get_col(v2,m,0);
gsl_matrix_set_row(m,1,v1);
gsl_matrix_set_col(m,1,v2);
KIMENET:
v1=[1 2 3]
v2=[1 4 7]
m=[[1 1 3] [1 4 3] [7 7 9]]

A gsl_matrix_get_[row/col](vektor,mátrix,i) függvény a vektorba másolja a mátrix i. [sorát/oszlopát]. A gsl_matrix_set_[row/col](mátrix,i,vektor) pedig a mátrix i. [sorába/oszlopába] másolja a vektort. Fontos, hogy megegyezzen a vektor mérete a sor vagy oszlopmérettel, különben hibajelzést kapunk.

KÓDRÉSZLET:
gsl_matrix_swap_rows(m,0,1);
KIMENET:
m=[[4 5 6] [1 2 3] [7 8 9]]

A gsl_matrix_swap_rows(mátrix,i,j)-vel felcserélhetjük a mátrix i. sorát és j. sorát. Ehhez hasonlóan a gsl_matrix_swap_cols(mátrix,i,j)-vel az oszlopokat cserélhetjük és a
gsl_matrix_swap_rowcol(mátrix,i,j)-vel pedig az i. sort a j. oszloppal.

KÓDRÉSZLET:
gsl_matrix_transpose(m);
KIMENET:
m=[[1 4 7] [2 5 8] [3 6 9]]

Értelemszerűen a gsl_matrix_transpose(mátrix) a mátrix transzponáltját állítja elő.

KÓDRÉSZLET:
gsl_matrix_view view;
gsl_matrix* v1;
view = gsl_matrix_submatrix(m,0,0,2,2);

v1 = &(view.matrix);
KIMENET:
v1=[[1 2] [4 5]]

A gsl_matrix_submatrix(mátrix,kezdő_sor,kező_oszlop,n,m) függvény a mátrixból a kezdő sor és oszloptól megadott nxm-es mátrixot vág ki. Ugyan az a helyzet mint a vektorok esetében, a kivágott mátrixot a megszokott módon átkonvertálhatjuk és használhatjuk, de a módosításra figyelni kell!

Hasonló a gsl_matrix_[row/col](mátrix,i) mely a mátrix i. [sorából/oszlopából] készít gsl_vector_view struktúrát melyet a megszokott módon kezelhetünk.

KÓDRÉSZLET:
gsl_vector_view view1 = gsl_matrix_diagonal(m);

gsl_vector_view view2 = gsl_matrix_subdiagonal(m,1);

gsl_vector_view view3 = gsl_matrix_superdiagonal(m,1);

gsl_vector* v1= &(view1.vector);

gsl_vector* v2= &(view2.vector);

gsl_vector* v3= &(view3.vector);
KIMENET:
v1=[1 5 9]
v2=[4 8]
v3=[2 6]

A mátrixból kivehetünk átlókat is. A gsl_matrix_diagonal(mátrix) a mátrix főátlóját, a gsl_matrix_subdiagonal(mátrix,i) a mátrix diagonálja alatti i. átlót, a
gsl_matrix_superdiagonal(mátrix,i) pedig a mátrix diagonálja fölötti i. átlót adja vissza vektorként.

Érdemes még megemlíteni a gsl_vector_view_array(tömb,n) függvényt, mely adott tömböt alakít gsl_vector_view struktúrára (n a vektor mérete) és a gsl_matrix_view_array(tömb,n,m) függvényt, mely a tömbből csinál egy nxm-es gsl_matrix_view struktúrájú mátrixot.

2.2. Műveletek

Legelőször nézzük meg a vektorműveleteket. Vektorok másolásához két függvény használatos, nevezetesen a gsl_vector_memcpy(cél_vektor,forrás_vektor) mely a cél_vektorba másolja a forrás_vektor-t, és a gsl_vector_swap(vektor1,vektor2) mely megcseréli a két vektor elemeit. Természetesen fontos, hogy a vektorok mindkét függvénynél megfelelő méretűek legyenek, különben hibajelzést kapunk.

Tekintsük át a két operandusú műveleteket. A következő függvényekkel lehet elvégezni az [összeadás/kivonás/szorzás/osztás] műveletét: gsl_vector_[add/sub/mul/div](vektor_a,vektor_b).

Fontos megjegyezni, hogy a művelet eredménye a legelső paraméterben megadott vektorban lesz eltárolva, felülírva így az ott lévő adatokat. Érdekes eredményt kapunk két vektor osztása esetén:

Legyen a=b=[1 1 1], c=[0 0 0], d=[5 5 5]

KÓDRÉSZLET:
gsl_vector_mul(a,c);
KIMENET:
a=[inf inf inf]
b=[inf inf inf]
c=[0 0 0]
d=[5 5 5]

Látható, hogy a GSL helyesen kezeli a nullával való osztást, illetve a végtelent pont ahogy elvárható és kipróbálható az is, hogy inf/inf=nan is igaz. Tehát a GSL ilyen esetekben nem ad vissza hibát, ezért érdemes odafigyelni rá. Azt, hogy egy vektor nullvektor e vagy sem azt a
gsl_vector_isnull(vektor) függvény segítségével ellenőrizhetjük. 1-et ad vissza ha nullvektor, nullát ha nem az.

Számos más használható művelet van még röviden ezek:

gsl_vector_scale(vektor,x) - a vektor elemeit szorozza x-el
gsl_vector_add_constant(vektor,x) - a vektor elemeihez hozzáad x-et
gsl_vector_max(vektor) - visszaadja a vektor legnagyobb elemét
gsl_vector_min(vektor) - visszaadja a vektor legkisebb elemét
gsl_vector_minmax(vektor,min,max) - min, max változókba teszi a maximális és minimális elemet
gsl_vector_max_index(vektor) - visszaadja a vektor legnagyobb elemének az indexét
gsl_vector_min_index(vektor) - visszaadja a vektor legkisebb elemének az indexét

gsl_vector_minmax_index(vektor,min,max) - min, max változókba teszi a maximális és minimális elem indexét

Térjünk át a mátrixműveletekre. Mátrixoknál is megtalálható a
gsl_matrix_memcpy(cél_mátrix,forrás_mátrix) és a gsl_matrix_swap(mátrix1,mátrix2). Szintén megtalálhatóak a gsl_matrix_[add/sub](mátrix_a,mátrix_b) [összeadás/kivonás] műveletek.
A [szorzás/osztás]-t a gsl_matrix_[mul/div]_elements(matrix_a,matrix_b) függvény végzi elemenként. Implementációra került a gsl_matrix_scale(mátrix,x) és a
gsl_matrix_add_constant(mátrix,x) is. És természetesen megtalálható az összes minimum és maximum keresést, nullmátrix ellenőrzést végző függvény mátrixos átirata is. Tehát összességében ugyanazon műveletek állnak rendelkezésre itt is mint a vektorok esetében és ugyanúgy hibajelzést kapunk nem megfelelő mátrix méretek esetén és ugyanúgy működik a nullmátrix-val való osztás is.

2.3. BLAS

Felmerülhet a kérdés, hogy hova lettek azok a műveletek, melyek két mátrix nem elemenkénti szorzását végzik, hol vannak azon függvények melyek mátrixok és vektorok szorzását valósítják mek, stb. Erre ad válasz a BLAS (Basic Linear Algebra Subprograms). Ez olyan függvények gyűjteménye melyek ezekre a problémákra vannak elkészítve, optimalizálva.

Három részre oszthatjuk a BLAS függvényeket:

1. Vektor műveletek:
2. Mátrix-vektor műveletek:
3. Mátrix-mátrix műveletek:

Legyen három mátrixunk:
A=[[0.11 0.12 0.13] [0.21 0.22 0.23]], B=[[1012 1022 1032] [1011 1021 1031]], C=[[0 0] [0 0]]

A következő műveletet szeretnénk elvégezni: C=A*B. Láthatjuk, hogy a szorzás elvégzéséhez az A vagy B mátrixot először transzponálni kellene, máskülönben nem lehet összeszorozni a két mátrixot egymással. Elárulhatjuk, hogy direkt választottuk így a két mátrixot. A megoldáshoz használhatjuk a gsl_blas_dgemm(trans_A,trans_B,alfa,mátrix_A,mátrix_B,béta,mátrix_C) függvényt. Ez a 3. típusú BLAS függvények közzé tartozik. A trans_[A/B] paraméter értéke a következőek lehetnek: CblasNoTrans, CblasTrans, CblasConjTrans. Ezzel végezzük el a megfelelő mátrixok transzponálását. Nézzünk meg egy példát:

KÓDRÉSZLET:
double a[] = ;
double b[] = ;
double c[] = ;

gsl_matrix_view A = gsl_matrix_view_array(a,2,3);
gsl_matrix_view B = gsl_matrix_view_array(b,2,3);
gsl_matrix_view C = gsl_matrix_view_array(c,2,2);

gsl_blas_dgemm(CblasNoTrans,CblasTrans,1.0,
&A.matrix,&B.matrix,0.0,&C.matrix);
KIMENET:
C=[[368.12 367.76] [674.72 674.06]]
B=[[1012 1022 1032] [1011 1021 1031]]

A példából látható, hogy a B mátrixot transzponáljuk, azonban az nem módosul véglegesen csak a számolás alatt lesz transzponált. Az eredmény természetesen a C mátrixban keletkezik. A képletben az alfát egyre kell állítani, a bétát pedig nullára, mivel nem akarjuk a C mátrixot belevenni a számolásba. Ezzel egyszerűen elvégezhetők a mátrixok szorzásai, s esetleg még más műveleteket is egybe tudunk fűzni egy függvényhívás alatt. Természetesen mint mindig, ha mégsem felelnének meg a mátrixok méretei a szorzás szabályának hibát kapunk a GSL-től. Megemlítjük még hogy az alfa és béta skalárszorzás elemenként értendő.

Nézzük meg mátrix vektorral való szorzását. Ez a 2. típusú BLAS függvények csoportjába tartozik, s a függvény a gsl_blas_dgemv(trans_A,alfa,mátrix_A,vektor_x,béta,vektor_y). Ennél a függvénynél is megadhatjuk az A mátrix transzponálását ugyanolyan módon, mint a fenti esetnél.

Most végezzük el két vektor szorzását: x=[1011 1012 1013] y=[0.11 0.12 0.13], keressük az x^Ty és az xy^T szorzatot. A várható eredmény első esetben egy szám, második esetben pedig egy mátrix.

KÓDRÉSZLET:
double x[] = ;
double y[] = ;
double p;
double Q[] =

gsl_vector_view v1 = gsl_vector_view_array(x,3);

gsl_vector_view v2 = gsl_vector_view_array(y,3);

gsl_blas_ddot(&v1.vector,&v2.vector,&p);
gsl_blas_dger(1.0,&v1.vector,&v2.vector,Q);
KIMENET:
p=367.76
Q=[[111.21 121.32 131.43] [112.31 122.52 132.73]
[113.41 123.72 134.03]]

A gsl_blas_ddot(vektor_x,vektor_y,eredmény) az x^Ty értékét számolja ki és azt az eredmény változóban adja vissza. A gsl_blas_dger(alfa,vektor_x,vektor_y,mátrix) a képlet alapján számolja ki az eredményt.

Ezen függvények mellett, természetesen számos más hasznos függvényt találhatunk a BLAS rutinok között, de azokat most nem részletezzük, minden pontosan megtalálható a GSL referencia dokumentációjában, mely letölthető a GSL weboldaláról.

2.4. Normák

Vektorokra vonatkozóan az euklideszi normát könnyedén meghatározhatjuk a
gsl_blas_snrm2(vektor) függvénnyel. Ez szintén a BLAS gyűjteményben található és a paraméterben megadott vektorhoz tartozó normát adja vissza visszatérési értékként. Sajnos az euklideszi normán kívül más norma számolásához nem biztosít függvényt a GSL.

Rosszabb a helyzet a mátrixok esetén. Egyáltalán nem találhatunk semmilyen beépített függvényt, mely mátrixokra normaszámítást tenne lehetővé. Ha valakinek mégis szüksége lenne rá, annak saját magának kell megírnia ezt a függvényt.

3. Transzformációk

3.1. QR felbontás

Vegyünk egy 3x4-es mátrixot és bontsuk fel egy ortogonális(Q) és egy felső trianguláris(R) mátrixra úgy hogy A=QR. Legyen A=[[3 0 2 -1] [0 -4 1 0] [4 0 5 3]].

KÓDRÉSZLET:
double a[] = ;

gsl_matrix_view A = gsl_matrix_view_array(a,3,4);
gsl_matrix* Q = gsl_matrix_alloc(3,3);
gsl_matrix* R = gsl_matrix_alloc(3,4);
gsl_vector* tau = gsl_vector_alloc(3);
gsl_matrix* C = gsl_matrix_alloc(3,4);

gsl_linalg_QR_decomp(&A.matrix,tau);
gsl_linalg_QR_unpack(&A.matrix,tau,Q,R);
gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,Q,R,0.0,C);
KIMENET:
A=[[3 0 2 -1] [0 -4 1 0] [4 0 5 3]]
Q=[[-0.6 0 -0.8] [0 1 0] [-0.8 0 0.6]]
R=[[-5 0 -5.2 -1.8] [0 -4 1 0] [0 0 1.4 2.6]]
C=[[3 0 2 -1] [0 -4 1 0] [4 0 5 3]]

A fenti kódrészlet tartalmazza a megoldást a felvetett kérdésre. A GSL-ben a QR felbontást két lépésben lehet elvégezni. Először az eredeti mátrixunkból a gsl_linalg_QR_decomp(mátrix,vektor) segítségével egy diagonális és felső trianguláris mátrixot és egy Householder vektort állítunk elő. Ezek után a gsl_linalg_QR_unpack(mátrix,vektor,mátrix_q,mátrix_r) függvénnyel az előbb kapott mátrixból és vektorból elő tudjuk állítani a szükséges Q illetve R mátrixokat. Fontos, hogy a Q mátrix mxm, míg az R mátrix mxn-es kell legyen. A tau Householder vektort pedig min(m,n) méretűre válasszuk. Amennyiben véletlenül nem ezt tennénk a számoláskor a GSL hibát fog jelezni, és kiírja nekünk, hogy melyik esetben nem stimmel a méret. Végül összeszoroztuk a QR mátrixokat és valóban az eredeti mátrixot kaptuk.

3.2. Householder transzformáció

Mint tudjuk Householder mátrix létrehozását a P=I-rvv^T képlettel tehetjük meg. A GSL közvetlenül a Householder mátrix előállítására nem tartalmaz függvényt, de a fentebb már említett függvények segítségével pillanatok alatt írhatunk egy kódot, mely mondjuk két vektor egymásba transzformálásához szükséges Householder mátrixot előállítja. Nézzük is ezt meg egy példával: Legyen v1=[1 2 3], v2=[3 2 -1]

KÓDRÉSZLET:
double a[] = ;
double b[] = ;
double dot;

gsl_vector_view v1 = gsl_vector_view_array(a,3);
gsl_vector_view v2 = gsl_vector_view_array(b,3);
gsl_vector_view v = gsl_vector_view_array(a,3);

gsl_matrix* H = gsl_matrix_alloc((&v1.vector)->size,
(&v1.vector)->size);
gsl_matrix* P = gsl_matrix_alloc((&v1.vector)->size,
(&v1.vector)->size);
gsl_matrix_set_zero(P);
gsl_matrix_set_identity(H);

gsl_vector_sub(&v.vector,&v2.vector);

gsl_blas_ddot(&v.vector,&v.vector,&dot);
gsl_blas_dger(2/dot,&v.vector,&v.vector,P);

gsl_matrix_sub(H,P);
gsl_blas_dgemv(CblasNoTrans,1.0,H,&v1.vector,0.0,&v.vector);
gsl_matrix_free(P);
gsl_matrix_free(H);
KIMENET:
H=[[0.6 0 0.8] [0 1 0] [0.8 0 -0.6]]
v=[3 2 -1]

A példában jól látható a megvalósítás módja. Legelőször beolvassuk a két vektort v1 és v2 vektorokba. Aztán kivonjuk őket egymásból, majd kiszámoljuk a vektor, ill. tenzoriális szorzatukat. A kapott eredményeket a Householder mátrix kiszámításához megadott képletbe behelyettesítve megkapjuk a H mátrixot. Ezután tesztképpen egy az első vektort megszorozzuk ezzel a mátrixxal és valóban eredményül a másik mátrixot kapjuk.

A fenti módszerrel könnyedén meghatározhatjuk más transzformálási esetekben is a Householder mátrixot.

4. Sajátérték, sajátvektor

A GSL lehetőséget ad mátrixoknál sajátértékek és sajátvektorok meghatározására is. Egyetlen kikötés van a mátrixra, hogy az szimmetrikus legyen. Nézzünk meg egy konkrét példát:

Legyen A=[[-1 2 3] [0 1 2] [-1 0 -1]].

KÓDRÉSZLET:
double a[] = ;

gsl_matrix_view A = gsl_matrix_view_array(a,3,3);
gsl_matrix* evec = gsl_matrix_alloc(3,3);
gsl_vector* eval = gsl_vector_alloc(3);

gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(3);

gsl_eigen_symmv(&A.matrix,eval,evec,w);

gsl_eigen_symmv_free(w);
KIMENET:
eval=[0 -2 1]
evec=[[0.7071 0.7071 0] [0 0 1] [-0.7071 0.7071 0]]

Minden számításnál kell egy workspace-t létrehozni. Ezt a gsl_eigen_symmv_alloc(n) függvénnyel tehetjük meg. Fontos, hogy csak nxn-es mátrixokkal dolgozhatunk, s az előbbi függvény n paramétere a mátrix méretével egyenlő. A sajátérték, sajátvektor meghatározásához a
gsl_eigen_symmv(mátrix,mátrix_sv,vektor_sj,workspace) függvényt használhatjuk. Az egész végeztével fel kell szabadítanunk a workspace-t a gsl_eigen_symmv_free(workspace)-vel.

Érdekességként megjegyezzük, hogy ha az eljárást nem szimmetrikus mátrixokra alkalmazzuk, semmilyen hibajelzést nem kapunk.

5. Interpolációk

Adva van néhány pont a egyazon síkban, s szeretnének egy olyan görbét kapni, melyet ezen pontok határoznak meg valamilyen módon. A keresett görbét interpolációval kaphatjuk meg. A GSL lehetőséget nyújt lineáris, polinominális és spline interpolációkra. Nézzünk egy konkrét példát:

KÓDRÉSZLET:
int i;
double xi,yi;
double x[10],y[10];

for (i = 0; i < 10; i++)


gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_interp *interp = gsl_interp_alloc(gsl_interp_polynomial,10);

gsl_interp_init(interp,x,y,10);
for(xi = x[0];xi < x[9];xi += 0.01)

gsl_interp_free(interp);
gsl_interp_accel_free(acc);

Kezdetben a kontrollpontokat az x és y tömbök tárolják. Itt most mi ezeket egy kezdőértékkel töltöttük fel. Eztán létrehozunk egy gsl_interp_accel struktúrát a gsl_interp_accel_alloc() függvénnyel, mely a tömbkezelést gyorsítja fel. Létrehozzuk magát az interpolációs struktúrát is (gsl_interp). A gsl_interp_alloc(interp_típus,méret) függvénnyel tehetjük ezt meg, melynek első paramétere mondja azt meg, hogy milyen interpolációt szeretnénk (gsl_interp_linear, gsl_interp_polinomial, gsl_interp_cspline). A paraméterlistában a méret az eredeti pontok számát jelenti. Ezek után végigmegyünk azon x tartományon, ahol ki szeretnénk számolni az interpolációs pontokat. A gsl_interp_eval(interpoláció,tömb_x,tömb_y,xi,accel) teszi ezt meg. Megadjuk neki az interpolációs struktúrát, az x és y koordinátákat tartalmazó tömböket, a keresett pont x koordinátáját (xi) és a gyorsító struktúrát. Visszatérési értékként kapjuk a megadott x-hez tartozó y-koordinátát amit felhasználhatunk a későbbiekben. A GSL nem tartalmaz semmilyen beépített ablakkezelést, hogy esetleg az interpoláció eredményét grafikusan lehessen ábrázolni, sajnos ezt saját magunknak kell megtenni.