走进BLAS/LAPACK（2）--blas

reference blas

任何事情都要讲究方式方法，只要思路和方法是对的，会有事半功倍之效，反之则会事倍功半。

我刚开始想要弄明白blas和lapack的使用方法，在百度和谷歌上搜索了很多文章，可是还是不得门而入。看到有人在测试函数中直接声明外接函数，也有人用头文件；有人使用fortran的接口，也有人用C的接口；写完了程序之后，用gcc手动链接，也有人用makefile，还有人调用静态库，等等。把人都弄晕了，不知从何下手，等我完全弄明白之后，回过头来看，发现它们本质是一样的。

blas是什么？归根结底只是用fortran写的一个程序包

（这里的blas单指reference blas，其它版本的implementation不一定是用fotran写的），那我想调用blas，就是想调用fortran的程序而已。因为个人不懂fortran，只了解一些简单的C/C++，所以对我来说，问题就归结于如何在C/C++中调用fortran程序，也就是C/C++与fortran的混合编程。想明白这一点，很多问题都迎刃而解。

安装

我使用的tarfile的版本号是3.6.0，下载的网址就是blas的官方网站netlib，cblas也可以在那里下载。

进入下载文件目录,后解压
$ tar -xzvf blas-3.6.0.tgz
会出现一个文件夹，里面包括了所有的fortran源文件和一个makefile文件和一个make.inc文件。其实到这里，已经可以用了，至于makefile和make.inc文件的作用等会再说。现在的问题就在于如何调用fortran源程序了

使用方法

ddot

举一个最简单的例子来说明：使用blas中的ddot函数来计算两向量相乘，d表示的是double类型，而dot表示的是点乘。如何来使用C调用fortran呢，这里有几点要注意。

• fortran中所有变量都表示的是指针，所以在C/C++中，输入变量都得是指针
• 在C/C++中声明调用的fortran函数时，在函数名称后要加一个下划线
• 在fortran中，矩阵是按照列来存储的，而在C/C++中矩阵是按行来存储的

第三点要注意一下，所谓按列或按行来存储，意味着在定义矩阵时，并不是动态生成二维数组，而是生成一维数组来存储矩阵，这一点很像python中处理矩阵的方式，用numpy来生成一维数组存储二维矩阵，不要搞错了。

下面是我写的一个调用ddot的C程序：dot.c

#include <stdio.h>

double ddot_(int *,double *,int *,double *,int *);
int main()
{
   int N=2,INCX=1,INCY=1;
   double X[2]={1.0,1.0};
   double Y[2]={2.0,2.0};
   double re;
   re=ddot_(&N,X,&INCX,Y,&INCY);
   printf("the result is:%f\n",re);
   return 0;
}

程序很简单，也可在声明ddot前加上extern表明它是一个外部函数，不过不加也不会出错。接下来就是C与fortran的混合编程了，这个就更简单了，我使用的是GNU套件，所以我用gcc来编译C程序，用gfortran编译fortran程序，最后用gcc来链接。

$ gcc -c dot.c -o dot.o
$ gfortran -c ddot.f -o ddot.o
$ gcc -o output dot.o ddot.o

为了方便，我把C程序和fortran程序都放在了同一个文件夹下面。执行一下二进制文件

$ ./output
the result is:4.000000

没有问题，结果是正确的。所有问题好像都结束了，不需要什么头文件，也不需要什么静态库，动态库，理论上这样做是没问题的，但实际上，这样做会累死人，为什么？看第二个例子

dgemm

在blas中有一套命名规则，dgemm中

• d表示的是double类型
• ge表示的是general，普通矩阵，不是对角，不是对称，也不是带状
• mm表示的是矩阵与矩阵相乘

更多的命名规则可以自行百度。gemm是一系列很有代表性的函数，因为它们经常被用来测试计算速度，来反映库的好坏。这里，我只用dgemm来测试，因为平常也主要是用double类型。

算法中的主要计算公式为
$C=\alpha AB+\beta C$

另外，还可以定义A和B是否为转置或共轭转置，输入参数一共有13个，详细的参数设置可以去官网上查询，这里主要介绍如何链接。C程序代码为：Cdgemm.c

#include <stdio.h>
#include <stdlib.h>
void dgemm_(char *,char *,int *,int *,int *,double *,double *,int *,double *,int *,double *,double *,int *);

int main()
{
    double a[4]={1.0,1.0,1.0,1.0},b[4]={2.0,2.0,2.0,2.0},c[4];
    char tra='n';
    char trb='n';
    int m=2,n=2,k=2,lda=2,ldb=2,ldc=2;
    double al=1.0,be=1.0;
    dgemm_(&tra,&trb,&m,&n,&k,&al,a,&lda,b,&ldb,&be,c,&ldc);
    printf("the matrix c is:%f,%f\n%f,%f\n",c[0],c[1],c[2],c[3]);

    return 0;
}

按照之前的方法先编译，再链接。

$ gcc -c Cdgemm.c -o Cdgemm.o
$ gfortran -c dgemm.f -o dgemm.o

在上述编译过程没有问题，开始链接

$ gcc -o output  dgemm.o Cdgemm.o

出错了，部分输出信息为：

dgemm.f:(.text+0xd3): undefined reference to `lsame_'
dgemm.f:(.text+0xf9): undefined reference to `lsame_'
dgemm.f:(.text+0x188): undefined reference to `lsame_'
dgemm.f:(.text+0x1b2): undefined reference to `lsame_'
dgemm.f:(.text+0x1f2): undefined reference to `lsame_'
dgemm.o:dgemm.f:(.text+0x21c): more undefined references to `lsame_' follow
dgemm.o: In function `dgemm_':
dgemm.f:(.text+0x2fb): undefined reference to `xerbla_'

原因是在dgemm.f中调用了lsame.f和xerbla.f这两个程序，而我没有编译，只能把lsame.f和xerbla,f复制到当前目录下，编译

$ gfortran -c lsame.f -o lsame.o
$ gfortran -c xerbla.f -o xerbla.o
$ gcc -o output dgemm.o Cdgemm.o lsame.o xerbla.o

又报错了

xerbla.o: In function `xerbla_':
xerbla.f:(.text+0x79): undefined reference to `_gfortran_st_write'
xerbla.f:(.text+0x90): undefined reference to `_gfortran_string_len_trim'
xerbla.f:(.text+0xb3): undefined reference to `_gfortran_transfer_character_write'
xerbla.f:(.text+0xd1): undefined reference to `_gfortran_transfer_integer_write'
xerbla.f:(.text+0xe0): undefined reference to `_gfortran_st_write_done'
xerbla.f:(.text+0xef): undefined reference to `_gfortran_stop_string'

解结办法有两个，一个用gfortran链接，另一个是在上面的链接最后加上 -lgfortran。

总结一下，要想得到最后的可执行文件，需要下面这些步骤:

$ gcc -c Cdgemm.c -o Cdgemm.o
$ gfortran -c dgemm.f -o dgemm.o
$ gfortran -c lsame.f -o lsame.o
$ gfortran -c xerbla.f -o xerbla.o
$ gfortran -o output dgemm.o Cdgemm.o lsame.o xerbla.o

当然，实际上是不需要写出这么多obj（目标文件）的，我写的这么细是为了说明哪一步会出问题。可以发现，自己手动链接十分麻烦，因为你不知道源代码中调用了什么其它函数，只能在链接报错的时候自己去找，这还只是一个简单的dgemm程序，如果是更复杂的程序，依赖关系很复杂，人为去链接根本不可能。

这个时候，库的作用就体现出来了，我只要把自己的代码与库链接在一起，编译器会自动把需要的obj文件链接在一起，根本不需要人为的关心。

上文中的make.inc和makefile在这里就要派上用场了。根据个人的配置修改好make.inc之后，make all一下就可以，会生成一个静态库，比如说叫libreblas.a，再执行下列命令

$ gcc -c Cdgemm.c -o Cdgemm.o
$ gfortran -o output Cdgemm.o ./libreblas.a

相比之下，简单了很多。