Пример использования процедуры xGEMM из библиотеки BLAS в программе на Си

BLAS (Basic Linear Algebra Subroutines) - библиотека базовых подпрограмм линейной алгебры. Это набор подпрограмм для выполнения базовых операций над векторами и матрицами. Библиотека BLAS имеет интерфейсы к языкам Fortran и C. Ниже приведен пример использования процедуры xGEMM из библиотека BLAS в программе на C.

Процедура xGEMM вычисляет следующее выражение: С = αAB + βC

где: A - матрица M*K

B - матрица K*N

C - матрица M*N

α, β - коэффициенты.

Первая буква в названии функции определяет тип данных элементов матрицы:

S - вещественный с одинарной точностью (например, SGEMM),

D - вещественный с двойной точностью (например, DGEMM),

С - комплексный с одинарной точностью (например, CGEMM),

Z - комплексный с двойной точностью (например, ZGEMM).

Прототип функции SGEMM в интерфейсе языка C выглядит следующим образом:

void cblas_sgemm(const enum CBLAS_ORDER Order,

const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB,

const int M, const int N,

const int K, const float alpha, const float *A,

const int lda, const float *B, const int ldb,

const float beta, float *C, const int ldc);

где: Order - определяет порядок следования элементов:

CblasRowMajor - матрицы хранятся по строкам (стандартно в C),

CblasColMajor - матрицы хранятся по столбцам;

TransA, TransB - определяет предварительные операции над матрицами A и B:

CblasNoTrans - ничего не делать,

CblasTrans - транспонировать,

CblasConjTrans - вычислить сопряженную матрицу;

M, N, K - размеры матриц;

alpha, beta - коэффициенты;

lda, ldb, ldc - число элементов в ведущей размерности матрицы (строке или столбце). Для массивов языка Си - число элементов в строке:

lda = K

ldb = N

ldc = N

Для использования процедур из библиотеки BLAS, необходимо, чтобы массивы в памяти хранились одним непрерывным блоком. Например, массив A из следующего примера нельзя передавать в качестве параметра в BLAS-процедуру:

float **A=(float**)malloc(M*sizeof(float*));

for (i=0;i<M;i++) A[i]=(float*)malloc(K*sizeof(float));

В нем не гарантировано, что строки будут следовать одна за другой. Следующий массив B использовать можно:

float **B=(float**)malloc(M*sizeof(float*));

float *B1=(float*)malloc(M*K*sizeof(float))

for (i=0;i<M;i++) B[i]=&B1[i*K];

Пример программы умножения двух матриц с использованием BLAS:

#include<stdio.h>

#include<cblas.h> // заголовочный файл C-интерфейса библиотеки BLAS

#define M 300

#define N 400

#define K 500

int main()

{ int i,j;

float A[M*K], // массив расположен в памяти одним непрерывным блоком

B[K][N], // этот тоже

*C; // и этот тоже будет непрерывным

C=(float*)malloc(M*N*sizeof(float)); // ну вот, непрерывный ;)

for (i=0;i<M;i++) // инициализируем массив A

for (j=0;j<K;j++)

A[i*K+j]=3*j+2*i;

for (i=0;i<K;i++) // инициализируем массив B

for (j=0;j<N;j++)

B[i][j]=5*i+j;

for (i=0;i<M*N;i++) C[i]=5; // инициализируем массив C (зачем?)

// перемножаем!!!

// нам нужно C=AB, поэтому сделаем beta=0

cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,

M,N,K,1.0,A,K,&B[0][0],N,0.0,C,N);

for (i=0;i<M;i++) // выводим результат на экран (наверно не влезет: 300*400)

{ for (j=0;j<N;j++) printf("%.2f ",C[i*N+j]);

printf("\n");

}

free(C);

return 0;

}

Библиотеку ATLAS (реализацию BLAS) для ccfit можно взять отсюда: atlas3.6.0_Linux_P4SSE2.tgz. Файл можно скачать программой wget и распаковать в домашний каталог (появится каталог Linux_P4SSE2).

При компиляции программы необходимо будет указать пути к каталогам include и lib и прилинковать необходимые библиотеки libcblas.a и libatlas.a

Пример команды компиляции:

gcc -I $HOME/Linux_P4SSE2/include -L $HOME/Linux_P4SSE2/lib -O3 prog.c -lcblas -latlas