BLAS i LAPACK

LAPACK je biblioteka rutina pisanih u FORTRANU 77 namijenjenih rješavanju linearnih algebarskih sustava, problema najmanjih kvadrata, računanju svojstvenih i singularnih vrijednosti matrice itd. Kod se može naći na adresi http://www.netlib.org/lapack/, a users' guide na adresi http://www.netlib.org/lapack/lug/index.html.

BLAS (Basic Linear Algebra Subprograms) je biblioteka rutina pisanih u FORTRANU 77 (adresa http://www.netlib.org/blas/) koje obavljaju osnovne algebarske operacije (kao npr. skalarni produkt) što efikasnije. Efikasnost tih rutina značajno ovisi o njihovoj usklađenosti s arhitekturom na kojoj se izvršavaju. Stoga je uputno koristiti ATLAS (Automatically Tuned Linear Algebra Software) koji automatski optimizira BLAS rutine na danoj platformi ( adresa http://math-atlas.sourceforge.net/).

Pozivanje LAPACKovih i BLASovih funkcija iz C-a i C++-a

Prilikom pozivanja LAPACKovih rutina iz C-a ili C++-a treba paziti na nekoliko stvari.

Program treba povezati (linkati) s odgovaraućim bibliotekama;
Sve rutine koje se pozivaju moraju biti ispravno deklarirane;
Argumenti moraju biti predani rutini onako kako to zahtijevaju programi pisani u FORTRANu 77.

Objasnimo poziv runtini pisanoj u FORTRANU 77 na primjeru rutine sgesv koja rješava realan sustav linearnih algebarskih jednadžbi.

Ako se pogleda u users' guide, http://www.netlib.org/lapack/lug/index.html, onda se može naći ovakav opis rutine:

SUBROUTINE SGESV( N, NRHS, A,	LDA, IPIV, B, LDB, INFO	)
      INTEGER	    INFO, LDA, LDB, N, NRHS
      INTEGER	    IPIV( * )
      REAL	    A( LDA, * ), B( LDB, * )
SVRHA
  SGESV	računa rješenje realnog sustava linearnih jednadžbi
     A * X = B,	gdje je A  N-x-N matrica a X i B su N-x-NRHS
  matrice.

  Koristi se LU dekompozicija s parcijalnim pivotiranjem  koja 
  faktorizira A u obliku
     A = P * L * U,
  gdje je P permutacijska matrca, L je donja trokutasta, a U je
  gornja trokutasta.  Faktorizirana forma se koristi za rješavanje
  sustava jednadžbi A *	X = B.

ARGUMENTI
  N	  (input) INTEGER
          Red matrice A,  N >= 0.
  NRHS	  (input) INTEGER
          Broj desnih strana, tj. redaka u matrici B.  NRHS >= 0.
  A	  (input/output) REAL polje, dimenzije (LDA,N)
          Na ulazu, N-x-N matrica A.  Na izlazu faktori L i U
          iz factorizacije A = P*L*U; jedinična dijagonala od faktora
          L se ne pamti.   
  LDA	  (input) INTEGER
          Vodeća dimenzija polja A. LDA >=max(1,N).
  IPIV	  (output) INTEGER polje, dimenzije (N)
          Pivotni indeksi koji definiraju permutacijsku matricu P; 
          redak i matrice je zamijenjen retkom IPIV(i).
  B	  (input/output) REAL polje, dimenzije (LDB,NRHS)
          Na ulazu, N-x-NRHS matrica desne strane B. Na izlazu, 
          ako je INFO= 0, N-x-NRHS matrica rješenja X.
  LDB	  (input) INTEGER
          Vodeća dimenzija polja  B. LDB >=max(1,N).
  INFO	  (output) INTEGER
           = 0:	uspješan izlaz 
           < 0:	INFO= -i, i-ti argument	ima nedozvoljenu vrijednost
           > 0:	INFO= i, U(i,i) je egzaktno nula. Faktorizacija	je 
           izračunata, ali je faktor U egzaktno singularan pa rješenje
           nije moglo biti izračunato.

Napomene

Imena kompiliranih fortranskih rutine obično su dana malim slovima, tako da bismo u C-u i C++-u pisali sgesv, a ne SGESV. Na nekim platformama (npr. linuxu) fortranski prevodilac imenu rutine dodaje underscore, dok C i C++ prevodioci to ne rade. Stoga tamo fortranski poziv
```
       call sgesv( .... )
```
prelazi u
```
       sgesv_( .... );
```
Fortranske tipove podataka treba zamijeniti odgovarajućim tipovima iz C-a (C++-a). Ekvivalencije su dane u sljedećoj tabeli:

REAL, REAL*4 float

REAL*8, DOUBLE PRECISION double

INTEGER int

COMPLEX, COMPLEX*8 float[2] (polje od 2 float-a)

COMPLEX*16 double[2] (polje od 2 double-a)
Fortran koristi prijenos svih argumenata po referenci. To znači da su svi formalni argumenti rutine pokazivači. Stoga ispred svakog skalarnog argumenta treba postaviti adresni operator. Na primjer, fortranski kod
```
       call sgesv(n, nrhs, A, lda, ipiv, ldb, info)
```
prelazi u
```
       sgesv_(&n, &nrhs, A, &lda, ipiv, B, &ldb, &info);
```
Neke LAPACKove rutine uzimaju argumente tipa string. U svim rutinama (osim testing i timing rutinama), samo prvi znak je važan. Stoga treba rutini dati pokazivač na znakovnu verijablu. Na primjer: fortranski poziv
```
   call dpotrf( 'Upper', n, a, lda, info )
```
prelazi u
```
   char s = 'U';
   dpotrf_(&s, &n, a, &lda, &info);
```
Između dvodimenzionalnih polja u FORTRANU i C, C++-u postoje velike razlike. Fortranska deklaracija
```
    REAL A(LDA, N)
```
definira polje od LDA*N float-ova u memoriji, spremljenih po stupcima: to znači da prvo idu elementi prvog stupca, zatim drugog itd. FORTRAN 77 (kao i C99) dozvoljava da se dimenzije dvodimenzionalnog polja rutini prenesu putem argumenata. Jedno dvodimenzionalno polje je samo pokazivač na prvi element polja. Stoga je dobar način alokacije memorije za matricu koja će biti predana LAPACKovoj rutini sljedeći:
```
    float *A;
    A = malloc( LDA*N*sizeof(float) ); 
```
Prilikom inicijalizacije matrice treba paziti na način na koji će fortranska rutina čitati matricu. Na primjer, kod koji bi inicijalizirao Hilbertovu matricu reda n
```
A =(a_ij),   a_ij = 1/(i+j-1)
```
izgledao bi ovako:
```
for(j=0; j<n; ++j){
        for(i=0; i<n; ++i){
            A[j*lda+i]=1.0/(i+j+1.0);  // A(i,j) 
        }
}
```
Indeks polja u FORTRANU ide od 1, dok u C/C++-u od 0, o čemu treba voditi računa.
Konačno, potrebno je napisati prototip rutine koju pozivamo. U našem slučaju to bi bilo
```
void sgesv_(int*, int*, float*, int*, int*, float*, int*, int*);
```
Sve fortranske SUBROUTINE treba deklarirati kao void. U C++-u deklaracija je
```
extern "C"{
void sgesv_(int*, int*, float*, int*, int*, float*, int*, int*);
}
```
Prilikom kompilacije programa potrebno je linkeru dati staze BLAS, LAPACK i fortranskih biblioteka. Ako su BLAS i LAPACK dobro instalirani bit će dovoljno navesti -lblas -llapack. Ime fortranske biblioteke može biti različito na različitim sustavima. Pod Linuxom treba zadati opciju -lf2c. Jedna naredba za kompilaciju pod Linuxom bi mogla biti
```
gcc ime_programa.c -lblas -llapack -lf2c -lm 
```
ili za C++ program
```
g++ ime_programa.cpp -lblas -llapack -lf2c -lm 
```
Poredak biblioteka je značajan.

Primjer poziva rutini sgesv na linux platformi iz jezika C dan je ovdje.

REAL, REAL*4	float
REAL*8, DOUBLE PRECISION	double
INTEGER	int
COMPLEX, COMPLEX*8	float[2] (polje od 2 float-a)
COMPLEX*16	double[2] (polje od 2 double-a)