图

BLAS 简介

预备知识 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 的相关页面

使用 cBLAS

预备知识 矩阵的储存

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

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

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

1
2
3
4
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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
// 计算矩阵—矢量乘法 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] << "  ";
    }
}

   编译与运行结果

1
2
3
4
5
6
7
8
$ g++ test_cblas.cpp -lblas
$ ./a.out
(1,2)  (2,3)  (3,4)

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

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

其他参数


1. 编程中的函数, 不是数学上的函数, 在一些编程语言(如 Fortran)中也叫子程序(subroutine)
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. 其实这个转置操作实现起来并不需要额外的运算, 在函数的代码中只需要把列主序(行主序)矩阵看作是行主序(列主序)的即可. 所以这么做比事先在内存中将矩阵元调换的方法快许多

致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条(需要权限) 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利