BLAS 简介

             

Prerequisite　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

Prerequisite　矩阵的储存

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

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

　　 类型

　　 级别（level）

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

　　 根据上述的类型和级别，BLAS 有一套系统的命名规则， 我们需要的函数名格式为 cblas_?gemv，其中 cblas_ 是固定前缀，问号表示 s, d, c, z, 中的一个，分别代表单精度（float），双精度（double），单精度复数和双精度复数4．ge 表示输入的矩阵是一个一般的矩阵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），m 和 n 分别指定矩阵的行数和列数，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)

其他参数

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 * 传递复数（任何指针都可以自动转换），所以只要内存中用两个连续的 float 或 double 表示一个复数即可
5. ^ 而不是对称矩阵或者厄米矩阵等，后者可以使用专门的函数使性能提升，例如 symv
6. ^ 其实这个转置操作实现起来并不需要额外的运算，在函数的代码中只需要把列主序（行主序）矩阵看作是行主序（列主序）的即可．所以这么做比事先在内存中将矩阵元调换的方法快许多

         

© 小时科技 保留一切权利