среда, 9 июня 2010 г.

Пример использования LAPACK, BLAS и ATLAS в Дебиан

Рассмотрим сначала использование LAPACK. Для начала создадим файл lapack.c следующего содержания:
$ cat lapack.c


/* Example C code for solving a linear system Ax=b using LAPACK */
#include <stdio.h>

/* Declare function prototype */
extern int sgesv_(int *n, int *nrhs, float *a, int *lda,
int *ipiv, float *b, int *ldb, int *info);

/* -- LAPACK driver routine (version 3.1) --
Purpose
=======

SGESV computes the solution to a real system of linear equations
A * X = B,
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

The LU decomposition with partial pivoting and row interchanges is
used to factor A as
A = P * L * U,
where P is a permutation matrix, L is unit lower triangular, and U is
upper triangular. The factored form of A is then used to solve the
system of equations A * X = B.

Arguments
=========

N (input) INTEGER
The number of linear equations, i.e., the order of the
matrix A. N >= 0.

NRHS (input) INTEGER
The number of right hand sides, i.e., the number of columns
of the matrix B. NRHS >= 0.

A (input/output) REAL array, dimension (LDA,N)
On entry, the N-by-N coefficient matrix A.
On exit, the factors L and U from the factorization
A = P*L*U; the unit diagonal elements of L are not stored.

LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,N).

IPIV (output) INTEGER array, dimension (N)
The pivot indices that define the permutation matrix P;
row i of the matrix was interchanged with row IPIV(i).

B (input/output) REAL array, dimension (LDB,NRHS)
On entry, the N-by-NRHS matrix of right hand side matrix B.
On exit, if INFO = 0, the N-by-NRHS solution matrix X.

LDB (input) INTEGER
The leading dimension of the array B. LDB >= max(1,N).

INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, U(i,i) is exactly zero. The factorization
has been completed, but the factor U is exactly
singular, so the solution could not be computed.

===================================================================== */

#define SIZE 3 /* dimension of matrix */

static int solve(float A[][SIZE], float *b) {

int i,j,info,n,nrhs,lda,ipiv[SIZE],ldb;
float AT[SIZE*SIZE];

/* Permute matrix */
for (i=0; i<SIZE; i++) {
for(j=0; j<SIZE; j++)
AT[j*SIZE+i]=A[i][j];
}

/* Invoke sgesv_ */
n = lda = ldb = SIZE; nrhs = 1;
sgesv_(&n, &nrhs, AT, &lda, ipiv, b, &ldb, &info);
return info;
}

int main(int argc, char **argv)
{
int i, j, pivot[SIZE], ok;
float A[SIZE][SIZE], b[SIZE];

/* Matrix A */
A[0][0]= 1.1; A[0][1]= 2.2; A[0][2]=-3.3;
A[1][0]= 4.4; A[1][1]=-5.5; A[1][2]= 6.6;
A[2][0]=-7.7; A[2][1]= 8.8; A[2][2]= 9.9;

/* Define right hand side vector */
b[0] = 0;
b[1] = 5.5;
b[2] = 11;

/* Call warpper function */
ok = solve(A, b);

/* Print out solution vector x */
for (j=0; j<SIZE; j++)
printf("%g\n", b[j]);
}



Компиллируем этот код с помощью команды:
$ gcc lapack.c -o lapack -llapack

Смотрим - какие пакеты нам нужны (просто для информации)
$ apt-file search /usr/lib/liblapack.so
lapack3: usr/lib/liblapack.so.3
lapack3: usr/lib/liblapack.so.3.0
lapack3-dev: usr/lib/liblapack.so


Рассмотрим теперь пример использования BLAS и еще один пример использования LAPACK в одном файле. Создадим файл blas.c такого содержания:
$ cat blas.c



#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <cblas.h>
#include <clapack.h>

int clapack_sgetrf(const enum CBLAS_ORDER Order, const int M, const int N, float *A, const int lda, int *ipiv);
// a is a column-major array of all the values in the matrix to invert
// The matrix's height and width are the same because it is a square matrix.
void invertMatrix(float *a, unsigned int height)
{
int info, ipiv[height];
info = clapack_sgetrf(CblasColMajor, height, height, a, height, ipiv);
info = clapack_sgetri(CblasColMajor, height, a, height, ipiv);
}

void displayMatrix(float *a, unsigned int height, unsigned int width)
{
int i, j;
for(i = 0; i < height; i++)
{
for(j = 0; j < width; j++)
{
printf("%1.3f ", a[height*j + i]);
}
printf("\n");
}
printf("\n");
}

void multiplyMatrix(float *a, unsigned int aheight, unsigned int awidth, float *b, unsigned int bwidth, float *c)
{
cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, aheight, bwidth, awidth, 1.0f, a, aheight, b, awidth, 0.0f, c, aheight);
}

int main(int argc, char *argv[])
{
int i;
float a[9], b[9], c[9];
srand(time(NULL));
for(i = 0; i < 9; i++)
{
a[i] = 1.0f*rand()/RAND_MAX;
b[i] = a[i];
}
displayMatrix(a, 3, 3);
invertMatrix(a, 3);
multiplyMatrix(a, 3, 3, b, 3, c);
displayMatrix(c, 3, 3);
return 0;
}

Компиллируем его:
$ gcc blas.c -o blasAndLapackUsage -llapack -L /usr/lib/atlas
Посмотрим дополнительную информацию про файлы, которые нужны для компиляции этого примера:
$ apt-file search /usr/include/clapack.h
atlas3-headers: usr/include/clapack.h
$ apt-file search /usr/include/cblas.h
refblas3-dev: usr/include/cblas.h
$ apt-file search /usr/lib/liblapack.so
lapack3: usr/lib/liblapack.so.3
lapack3: usr/lib/liblapack.so.3.0
lapack3-dev: usr/lib/liblapack.so

$ apt-file search /usr/lib/atlas
atlas3-base: usr/lib/atlas/libblas.so.3
atlas3-base: usr/lib/atlas/libblas.so.3.0
atlas3-base: usr/lib/atlas/liblapack.so.3
atlas3-base: usr/lib/atlas/liblapack.so.3.0
atlas3-base-dev: usr/lib/atlas/libblas.a
atlas3-base-dev: usr/lib/atlas/libblas.so
atlas3-base-dev: usr/lib/atlas/liblapack.a
atlas3-base-dev: usr/lib/atlas/liblapack.so

Другие полезные ссылки по теме:
ATLAS FAQ
Example C code for solving a linear system Ax=b using LAPACK
Compiling Programs Containing LAPACK Routines
LAPACK build and test guide
Using Lapack Routines in C Programs
netlib, Frequently Asked Questions
BLAS Frequently Asked Questions (FAQ)
GNU Scientific Library -- Reference Manual
Driver Routines
The Linear Algebra Package

Комментариев нет: