走进BLAS/LAPACK(2)--blas

时间:2022-07-17 04:08:05

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=αAB+β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

相比之下,简单了很多。