BLAS 简介

                     

贡献者: addis; junyi

预备知识 C 语言基础

   BLAS(basic linear algebra subroutine) 是一系列基本线性代数运算函数1接口(interface)标准。 这里的线性代数运算是指例如矢量的线性组合,矩阵乘以矢量,矩阵乘以矩阵等。接口在这里指的是诸如哪个函数名实现什么功能,有几个输入和输出变量,分别是什么。

   BLAS 被广泛用于科学计算和工业界,已成为业界标准。在更高级的语言和库中,即使我们不直接使用 BLAS 接口,它们也是通过调用 BLAS 来实现的(如 Matlab 中的各种矩阵运算)。

   BLAS 原本是用 Fortran 语言写的,但后来也产生了 C 语言的版本 CBLAS,接口与 Fortran 的略有不同(例如使用指针传递数组),但大同小异。

   注意 BLAS 是一个接口的标准而不是某种具体实现(implementation)。简单来说,就是不同的作者可以各自写出不同版本的 BLAS 库,实现同样的接口和功能,但每个函数内部的算法可以不同。 这些不同导致了不同版本的 BLAS 在不同机器上运行的速度也不同。

   BLAS 的官网是 Netlib,可以浏览完整的说明文档以及下载源代码。这个版本的 BLAS 被称为 reference BLAS,运行速度较慢,通常被其他版本用于衡量性能。对于 Intel CPU 的计算机,性能最高的是 Intel 的 MKL(Math Kernel Library) 中提供的 BLAS。MKL 虽然不是一个开源软件,但目前可以免费下载使用。如果想要免费开源的版本,可以尝试 OpenBlas 或者 ATLAS2。另外,无论是否使用 MKL,BLAS 的文档都推荐看 MKL 的相关页面

1. 使用 CBLAS

预备知识 矩阵的储存

   我们这里介绍一个 CBLAS 的例子3。我们来测试矩阵和矢量的乘法函数 gemv,参考文档见这里

   BLAS 可以分为四种类型和三个级别

   类型

   级别(level)

  1. 向量与向量操作(复杂度 $n$)
  2. 矩阵与向量操作(复杂度 $n^2$)
  3. 矩阵与矩阵操作(复杂度 $n^3$)

   根据上述的类型和级别,BLAS 有一套系统的命名规则, 我们需要的函数名格式为 cblas_?gemv,其中 cblas_ 是固定前缀,问号表示 s, d, c, z, 中的一个,分别代表单精度(float),双精度(double),单精度复数和双精度复数4ge 表示输入的矩阵是一个一般的矩阵5,以行主序或者列主序线性储存,mv 表示矩阵与向量操作。

   我们接下来以 cblas_zgemv 为例,先来看函数声明。

void cblas_zgemv (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE trans,
const MKL_INT m, const MKL_INT n, const void *alpha, const void *a,
const MKL_INT lda, const void *x, const MKL_INT incx, const void *beta,
void *y, const MKL_INT incy);

   BLAS 接口给人的第一感觉就是冗长,为什么实现一个简单的功能需要这么多变量?因为这个接口具有相当大的灵活性。例如可以使用一个列主序(行主序)矩阵的子矩阵作为矩阵,又例如可以使用一个列主序(行主序)矩阵的某一行(列)作为矢量;例如可以在做乘法以前对矩阵进行转置6;又例如可以把相乘的结果累加到输出矢量已有的值上,而不是直接覆盖。

   作为一个最简单的例子,我们可以用以下测试函数,在这个例子中,我们只需要关心几个变量,即 Layout(一个 enum 类型)指定行主序(CblasRowMajor)还是列主序(CblasColMajor),mn 分别指定矩阵的行数和列数,a, x, y 分别是矩阵第一个元的指针。

代码 1:blas_demo.cpp
// 计算矩阵—矢量乘法 y = a * x
#include <iostream>
#include <complex>
#include <cblas.h>
int main() {
    using namespace std;
    typedef complex<double> Comp; // 定义复数类型
    int Nr = 2, Nc = 3; // 矩阵行数和列数
    // 矩阵和矢量分配内存
    Comp *a = new Comp [Nr*Nc];
    Comp *x = new Comp [Nc];
    Comp *y = new Comp [Nr];
    Comp alpha(1, 0), beta(0, 0);
    // x 矢量赋值
    for (int i = 0; i < Nc; ++i) {
        x[i] = Comp(i+1., i+2.);
    }
    // a 矩阵赋值
    for (int i = 0; i < Nr*Nc; ++i) {
        a[i] = Comp(i+1., i+2.);
    }
    // 做乘法
    cblas_zgemv(CblasColMajor, CblasNoTrans, Nr, Nc, &alpha, a,
        Nr, x, 1, &beta, y, 1);
        
    // 命令行分别输出 x, a, y
    for (int i = 0; i < Nc; ++i) {
        cout << x[i] << "  ";
    }
    cout << "\n" << endl;
    for (int i = 0; i < Nr; ++i) {
        for (int j = 0; j < Nc; ++j) {
            cout << a[i + Nr*j] << "  ";
        }
        cout << endl;
    }
    cout << "\n" << endl;
    for (int i = 0; i < Nr; ++i) {
        cout << y[i] << "  ";
    }
}

   编译与运行结果

$ g++ test_cblas.cpp -cblas
$ ./a.out
(1,2)  (2,3)  (3,4)

(1,2)  (3,4)  (5,6)
(2,3)  (4,5)  (6,7)

(-18,59)  (-21,74)

64 位整数

   注意使用 64bit 时,头文件不需要改变,而是在 #include <cblas.h> 之前定义 #define CBLAS_INT long long

   另外在 ubuntu 上需要安装 libblas64-dev

其他参数

2. XBLAS

3. 从源码编译 CBLAS

   一般用 (C)BLAS 也需要用 LAPACK(E),所以可以直接用 LAPACK 的源码编译,详见 “Lapack 笔记”。

4. OpenBlas


1. ^ 编程中的函数,不是数学上的函数,在一些编程语言(如 Fortran)中也叫子程序(subroutine)。可以简单地认为 Fortran 函数有返回值,子程序没有,如下面的例子中,传进和传出的是指针
2. ^ 至于安装,在 Windows 系统上作者推荐在 Visual Studio 的基础上安装 MKL 或者 Parallel Studio,这些软件都比较大,可能需要较长时间下载安装。在 Linux 系统上,可以直接用 apt-get 等安装开源版本,也可以下载 intel MKL 的 deb 安装包按照提示安装。
3. ^ 要在 Linux 上编译,见 “在 Linux 上编译第一个 C++ 程序
4. ^ 注意 C 语言中内置的复数类型和 C++ 标准库中的 std::complex<> 不同,但由于 CBLAS 用指针 void * 传递复数(任何指针都可以自动转换),所以只要内存中用两个连续的 floatdouble 表示一个复数即可
5. ^ 而不是对称矩阵或者厄米矩阵等,后者可以使用专门的函数使性能提升,例如 symv
6. ^ 其实这个转置操作实现起来并不需要额外的运算,在函数的代码中只需要把列主序(行主序)矩阵看作是行主序(列主序)的即可。所以这么做比事先在内存中将矩阵元调换的方法快许多

                     

© 小时科技 保留一切权利