CBLAS的安装与使用

CBLAS是BLAS的C语言接口。BLAS的全称是Basic
Linear Algebra Subprograms，中文大概可以叫做基础线性代数子程序。主要是用于向量和矩阵计算的高性能数学库。本身BLAS是用Fortran写的，为了方便C/C++程序的使用，就有了BLAS的C接口库CBLAS。BLAS的主页是http://www.netlib.org/blas/，CBLAS的下载地址也可以在这个页面上找到。

 

CBLAS安装需要先装BLAS，从主页上下载blas.tgz，解压，根据系统修改make.inc和Makefile，make，就会生成一个blas_LINUX.a文件。然后，下载cblas.tgz，解压，在目录下将Makefile.*文件改名或者做一个链接文件为Makefile.in文件，比如在linux下就是ln
-s Makefile.LINUX Makefile.in，根据具体情况修改Makefile.in文件，主要是BLAS的库文件路径BLLIB和CBLAS的安装目录CBDIR，make
help就可以打印出可以使用的make命令，要生成全部文件就是用make all。在$(CBDIR)目录下的$(CBLIBDIR)将生成CBLAS的库文件$(CBLIB)，cblas_LINUX.a。

 

在CBLAS的安装目录$(CBDIR)下的src目录中有个cblas.h是包括的CBLAS的函数和常量的头文件，使用CBLAS的时候就需要这个头文件，同时还需要BLAS的库文件$(BLLIB
)和CBLAS的库文件$(CBLIB)。

 

CBLAS/BLAS分为3个level，level1是用于向量的计算，level2是用于向量和矩阵之间的计算，level3是矩阵之间的计算。比如计算矩阵的乘法就是属于level3，这里就用矩阵乘法来学习使用CBLAS。

 

计算矩阵乘法的函数之一是 cblas_sgemm，使用单精度实数，另外还有对应双精度实数，单精度复数和双精度复数的函数。在此以 cblas_sgemm为例。

 

函数定义为：

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)

关于此函数的详细定义可以在http://www.netlib.org/blas/sgemm.f找到，只不过是fortran语言的,这个C语言版的略有差别。

此函数计算的是 C = alpha*op( A )*op( B ) + beta*C,

 

const enum CBLAS_ORDER Order，这是指的数据的存储形式，在CBLAS的函数中无论一维还是二维数据都是用一维数组存储，这就要涉及是行主序还是列主序，在C语言中数组是用行主序，fortran中是列主序。我还是习惯于是用行主序，所以这个参数是用CblasRowMajor，如果是列主序的话就是CblasColMajor。

 

const enum CBLAS_TRANSPOSE TransA和                 const enum CBLAS_TRANSPOSE TransB，这两个参数影响的是op(
A )和op( B)，可选参数为CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113，其中TransA
= CblasNoTrans, op( A ) = A，TransA = CblasTrans, op( A ) = A'，TransA = CblasConjTrans, op( A ) = A'。 TransB类似。

 

const int M，矩阵A的行，矩阵C的行

const int N，矩阵B的列，矩阵C的列

const int K，矩阵A的列，矩阵B的行

const float alpha， const float beta，计算公式中的两个参数值，如果只是计算C=A*B，则alpha=1,beta=0

 

const float *A， const float *B， const float *C，矩阵ABC的数据

 

const int lda， const int ldb， const int ldc，在BLAS的文档里，这三个参数分别为ABC的行数，但是实际使用发现，在CBLAS里应该是列数。

 

我在这里计算两个简单矩阵的乘法。

A:

1,2,3

4,5,6

7,8,9

8,7,6

 

B:

5,4

3,2

1,0

 

程序代码：

//因为程序是C++，而CBLAS是C语言写的，所以在此处用extern关键字

extern"C"

{

    #include<cblas.h>  

}

#include<iostream>

using namespace std;

int main(void)
{

   

    const enum CBLAS_ORDER Order=CblasRowMajor;

    const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;

    const enum CBLAS_TRANSPOSE TransB=CblasNoTrans;

    const int M=4;//A的行数，C的行数

    const int N=2;//B的列数，C的列数

    const int K=3;//A的列数，B的行数

    const float alpha=1;

    const float beta=0;

    const int lda=K;//A的列

    const int ldb=N;//B的列

    const int ldc=N;//C的列

    const float A[K*M]={1,2,3,4,5,6,7,8,9,8,7,6};

    const float B[K*N]={5,4,3,2,1,0};

    float C[M*N];

   

    cblas_sgemm(Order, TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);

     

    for(int i=0;i<M;i++)

    {

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

       {

           cout<<C[i*N+j]<<"/t";

       }

       cout<<endl;

    }

   

    return EXIT_SUCCESS;

}

在编译的时候需要带上cblas_LINUX.a和blas_LINUX.a，比如，

g++ main.cpp cblas_LINUX.a blas_LINUX.a -o main

当然，这里假定是这两个.a文件是放在可以直接访问的位置，或者写全路径也可以。

 

这种做法在CentOS.5下顺利通过，但是在我的Ubuntu.7.10下出了问题，blas_LINUX.a正常编译生成，但在链接的时候出了错误，所以只好从源里安装了atlas，sudo
apt-get install atlas3-base，在/usr/lib/atlas/目录下就会有libblas.*和liblapack.*库文件，只需要在链接的时候用这里的blas库文件替换上文安装的BLAS就可以正常编译通过。

 

另外，在GSL下也有BLAS和CBLAS，在boost里有ublas也提供CBLAS/BLAS的功能，有时间也拿来研究研究。

另附GNU上的一个例子：

The following program computes the product of two matrices using the Level-3 BLAS function SGEMM,

[ 0.11 0.12 0.13 ] [ 1011 1012 ] [ 367.76 368.12 ]
[ 0.21 0.22 0.23 ] [ 1021 1022 ] = [ 674.06 674.72 ]
[ 1031 1032 ]
The matrices are stored in row major order but could be stored in column major order if the first argument of the call to cblas_sgemm was changed to CblasColMajor.

#include <stdio.h>
#include <gsl/gsl_cblas.h>

int
main (void)
{
int lda = 3;

float A[] = { 0.11, 0.12, 0.13,
0.21, 0.22, 0.23 };

int ldb = 2;

float B[] = { 1011, 1012,
1021, 1022,
1031, 1032 };

int ldc = 2;

float C[] = { 0.00, 0.00,
0.00, 0.00 };

/* Compute C = A B */

cblas_sgemm (CblasRowMajor,
CblasNoTrans, CblasNoTrans, 2, 2, 3,
1.0, A, lda, B, ldb, 0.0, C, ldc);

printf ("[ %g, %g\n", C[0], C[1]);
printf (" %g, %g ]\n", C[2], C[3]);

return 0;
}
To compile the program use the following command line,

$ gcc -Wall demo.c -lgslcblas
There is no need to link with the main library -lgsl in this case as the CBLAS library is an independent unit. Here is the output from the program,

$ ./a.out
[ 367.76, 368.12
674.06, 674.72 ]