「C++」 线性运算库汇总
1 引言
MKL 和 OpenBLAS 都提供了比较好的性能，MKL 性能还更好一点；
2 常用库1
2.1 EIGEN
是一个线性算术的开源 C++ 模板库；功能强大、快速、优雅以及支持多平台，可以使用该库来方便处理一些矩阵的操作，达到类似 matlab 那样的快捷；
支持稠密矩阵(Dense Matrix)和稀疏矩阵(Sparse Matrix)运算，并内置实现了对两种矩阵的求解器(Solver)；使用 Eigen 十分简单，因为是模板库，所以只需要将头文件拷入 include 目录即可，缺点是编译时比较费时。Eigen 的稠密矩阵运算在某些平台上的性能甚至接近 Intel MKL；虽然 Eigen 的求解器从效率来讲并不算最好，但将其作为矩阵容器是一个很好的选择，在求解大型问题时再通过连向诸如 SuiteSparse 等顶级求解器即可；TensorFlow 即使用了Eigen；在 3.3 之后，Eigen 引入了对其他 BLAS 的接口，支持了包括 MKL，OpenBLAS 等高性能库，使得其性能还可能更进一步提高；
需要定义 EIGEN_NO_DEBUG 阻止运行时 assertion；编译单线程版本需要开启 -DEIGEN_DONT_PARALLELIZE. 在试验中，我们采用 EIGEN 原生 BLAS 实现；
2.2 Intel MKL
英特尔数学核心函数库是一套经过高度优化和广泛线程化的数学例程，专为需要极致性能的科学、工程及金融等领域的应用而设计；它可以为当前及下一代英特尔处理器提供性能优化，包括更出色地与 Microsoft Visual Studio、Eclipse和XCode相集成；英特尔 MKL 支持完全集成英特尔兼容性 OpenMP 运行时库，以实现更出色的 Windows/Linux 跨平台兼容性；在试验中的多线程版本需要链接到 mkl_gnu_thread，而不是 mkl_intel_thread，单线程版本需要链接到 mkl_sequential_thread；
2.3 OpenBLAS
是一个高性能多核 BLAS 库，是 GotoBLAS2 1.13 BSD 版本的衍生版；OpenBLAS 的编译依赖系统环境，并且没有原生单线程版本，在实验这哦那个，通过设置 OMP_NUM_THREADS=1 来模拟单线程版本，可能会带来一点点的性能下降；
2.4 Armadillo
Armadillo 是一款年轻的线性运算模板库，使用了模板元编程技术，拥有友好的交互接口，十分贴近 matlab 风格，迭代速度也很快；并且可以使用其他 BLAS 实现如 OpenBLAS、MKL 和 AND ACM（ACML） 作为计算后端，实现稠密矩阵运算的加速；但是 Armadillo 的稀疏矩阵求解只能连向 SuperLU，不支持连接 SuiteSparse（Cholmod 性能更好）；
3 稀疏矩阵求解库
Eigen 与 Armadillo 在对于一些特定问题的求解并不是最优，比如高性能计算界研究了很多年的稀疏矩阵求解(Sparse Matrix Solver)；对于标准的 LU/QR/Cholesky 分解，高性能计算界已经有成熟高效的计算库如 SuiteSparse；所以在追求计算速度的科研领域，还是需要掌握这些更复杂的求解库的使用；
3.1 SuiteSparse
SuiteSparse 是现在科学界和工业界最强大的基于 BLAS/LAPACK 的多种稀疏矩阵求解算法合集，其中包括了LU分解算法(UMFPACK)， Cholesky分解(CHOLMOD)， QR分解算法(SPQR)等算法的实现；Matlab 的矩阵求解就是使用的 SuiteSparse 的实现，可以说是相当的强大；笔者工作中就极度依赖 SuiteSparse 提供的加速，比如使用 Eigen 的自带稀疏矩阵求解器就要比使用 SuiteSparse 慢一两个数量级；
3.2 SuperLU
SuperLU 针对 LU 分解所写的库，和 UMFPACK 具有相同的使用范围，在各测试集上表现也各有优劣，笔者较为推荐 UMFPACK，因为可以在 SuiteSparse 中一并安装；
4 术语
在涉及线性运算的领域，几个术语 BLAS, LAPACK 是绕不过的；
4.1 BLAS LAPACK 的含义与历史
BLAS 全称为 Basic Linear Algebra Subprograms，是一个API标准，其设计目的是为了提供对向量和矩阵计算的基础支持，具体分为BLAS1(向量间运算)，BLAS2（矩阵与向量计算），BLAS3（矩阵与矩阵计算）；设计 BLAS 的原因是由于不同的计算机硬件架构（总线、各级缓存、CPU寄存器数量），逻辑上相同的运算在写为机器码过程中有很大的优化空间，所以通过这样一个标准来分离底层运算及更高级的调用，达到较好的跨平台性及运算速度；现在较为知名的 BLAS 实现是 Netlib 的 cBLAS， 第一个做到跨平台可伸缩的 ATLAS，boost 提供的 ublas，以及 GotoBLAS2 的继承改进版 OpenBLAS，还有Intel专门为Intel处理器开发的MKL；
LAPACK 全称 Linear Algebra Package，是一个基于 Level 3 BLAS 的线性代数数值计算库，其继承了 LINPACK 的线性方程求解和最小二乘求解，EISPACK的特征值求解，以及实现了多种矩阵分解算法（LU, Cholesky, QR, SVD, Schur, generalized Schur）。最新的 LAPACK 是由Fortran 90写成的。LAPACKE 是 LAPACK 的 C 语言API；
4.2 现有的BLAS库实现
比较全的介绍可以在 BLAS 的 wiki 中找到，以下主要谈一谈笔者的体会；
OpenBLAS 是笔者现在使用的 BLAS 库，其继承了 GotoBLAS2 的设计，还对 LAPACK 实现进行了优化，具有较好的性能。虽然赶不上 MKL 在 Intel X86-64 计算机上的性能，但是由于其开源性，笔者可以不用担心哪天 MKL 又不免费了。而且正是因为这些开源软件不断提高的品质才使得商业公司不得不免费提供他们的产品给社区以保持市场占有。所以虽然 OpenBLAS 在性能上落于 MKL 之下，笔者仍然将其放在第一位；
MKL 在其主页打出了”the fastest”的标语，的确是底气很足；由于是 Intel 的亲儿子，MKL 获得了很多内部支持，能够优化到CPU的最底层；在 Intel CPU 上没有比它更快的 BLAS 实现了，而且好消息是现在 MKL 可以免费下载用于学习科研目的；
ATLAS 全称 Automatically Tuned Linear Algebra Software，是一个可以对任何计算机硬件生成优化过的BLAS库的软件；虽然在对指定硬件架构上性能不及专门优化过的BLAS库，其自动生成能力使得它能够对任意新系统生成第一个优化过的BLAS库，并作为后来者的对比基准；
cBLAS 是 Netlib 提供的基于 Fortrain 实现的 BLAS 的 c 语言API，其存在价值在于通常默认安装于系统中；
uBLAS 是 boost 包提供的 C 语言 BLAS 实现，性能不敢恭维；
cuBLAS 是 Nvidia 为其 GPU 提供的线性运算库，包含于 CUDA 中，也可以单独下载安装；由于稠密矩阵运算是高度并行的任务，所以 GPU 更擅长做基础向量运算；在 Nvidia 官方给的测试中，使用 cuBLAS 的 Tesla P100 比使用 MKL 的 Xeon Dual E5-2699 v4 快了5倍；
clBLAS 是 BLAS 基于 OpenCL 的实现，主要针对并行和 AMD GPU进行了优化；clBLAS的优点是它为开源项目，其作用之于 Nvidia 类似于 OpenBLAS 之于 Intel；
5 总结
5.1 选型
- 对比
备选：MKL、OpenBLAS、Eigen、Armadillo；
接口易用程度：Eigen > Armadillo > MKL >/OpenBLAS;
速度：MKL = OpenBLAS > Eigen(with MKL) > Eigen > Armadill;
- OpenBLAS 没有单核版本；
- MKL 学习成本高，但性能最强；
- Armadillo 效率和接口易用性不如 Eigen；
- Eigen 原生实现速度不如 MKL、OpenBLAS，但是使用 MKL 作为后台，性能和 MKL 原声实现几乎一样；
附录
A 参考资料
B 性能对比2
- 每个测试程序的编译都采用 “-O4 -msse2 -msse3 -msse4” 优化， 通过设置 OMP_NUM_THREADS 来控制使用的线程数量. 除了 OpenBLAS，其他两个库的测试程序都分别有单线程和多线程的编译版本；*
单核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 500 | 0.04159 | 0.03122 | 0.03058 |
| 1000 | 0.31789 | 0.24339 | 0.23730 |
| 1500 | 1.04589 | 0.81445 | 0.79869 |
| 2000 | 2.37567 | 1.92036 | 1.87102 |
| 2500 | 4.68266 | 3.78569 | 3.64548 |
| 3000 | 8.28073 | 6.42630 | 6.29797 |
| 3500 | 13.07470 | 10.25096 | 9.98417 |
| 4000 | 19.34550 | 15.21931 | 14.87500 |
| 4500 | 27.52767 | 21.45024 | 21.18227 |
| 5000 | 37.67552 | 29.31631 | 29.07229 |
8 核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 1000 | 0.05658 | 0.03955 | 0.06468 |
| 2000 | 0.34981 | 0.26200 | 0.23879 |
| 3000 | 1.20781 | 0.85449 | 0.80737 |
| 4000 | 2.65490 | 1.90273 | 1.88366 |
| 5000 | 5.03304 | 3.73005 | 3.67966 |
| 6000 | 8.78654 | 6.52766 | 6.31980 |
| 7000 | 13.55611 | 10.13758 | 10.07120 |
| 8000 | 19.81634 | 15.03530 | 14.89440 |
| 9000 | 29.11329 | 21.54359 | 21.26992 |
| 10000 | 39.01563 | 29.93075 | 29.22034 |
16 核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 1000 | 0.05708 | 0.02185 | 0.03897 |
| 2000 | 0.26694 | 0.13807 | 0.30461 |
| 3000 | 0.70686 | 0.43692 | 0.93511 |
| 4000 | 1.45129 | 0.97720 | 2.06761 |
| 5000 | 2.59477 | 1.90665 | 2.49280 |
| 6000 | 5.43438 | 3.30945 | 7.01299 |
| 7000 | 8.01124 | 5.17896 | 6.84496 |
| 8000 | 11.22280 | 7.81439 | 12.99240 |
| 9000 | 15.15625 | 11.08906 | 21.82488 |
| 10000 | 19.91151 | 15.22039 | 30.86908 |
32 核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 1000 | 0.04003 | 0.02792 | 0.02244 |
| 2000 | 0.51213 | 0.14363 | 0.16990 |
| 3000 | 1.13647 | 0.51105 | 0.54635 |
| 4000 | 1.58793 | 1.10219 | 1.26401 |
| 5000 | 2.88341 | 2.07923 | 2.48735 |
| 6000 | 5.92779 | 3.42785 | 4.26794 |
| 7000 | 7.91650 | 5.32176 | 6.69391 |
| 8000 | 11.96467 | 7.65395 | 9.98951 |
| 9000 | 17.45420 | 10.28328 | 14.14108 |
| 10000 | 23.31314 | 15.10077 | 19.34171 |
40 核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 1000 | 0.03691 | 0.02877 | 0.01779 |
| 2000 | 0.37739 | 0.14037 | 0.13655 |
| 3000 | 0.61183 | 0.41057 | 0.44113 |
| 4000 | 2.43670 | 1.02625 | 1.01414 |
| 5000 | 3.18099 | 1.91092 | 1.97898 |
| 6000 | 8.24002 | 2.96157 | 3.40685 |
| 7000 | 11.59889 | 4.68312 | 5.38634 |
| 8000 | 9.50613 | 6.98434 | 7.95971 |
| 9000 | 14.83066 | 9.60891 | 11.37585 |
| 10000 | 23.67187 | 15.52151 | 15.52680 |
48 核
| Matrix-Dimension | Eigen | MKL | OpenBLAS |
|---|---|---|---|
| 1000 | 0.03635 | 0.02398 | 0.01548 |
| 2000 | 0.36417 | 0.13408 | 0.11496 |
| 3000 | 2.32388 | 0.39291 | 0.36669 |
| 4000 | 2.32030 | 1.13244 | 0.85790 |
| 5000 | 2.08269 | 1.75812 | 1.66785 |
| 6000 | 8.70766 | 2.98694 | 2.85609 |
| 7000 | 8.23543 | 4.62340 | 4.53257 |
| 8000 | 21.18603 | 6.68886 | 6.72820 |
| 9000 | 19.86504 | 9.59635 | 9.50597 |
| 10000 | 16.10920 | 13.13038 | 13.04432 |
-
. 矩阵运算库blas, cblas, openblas, atlas, lapack, mkl之间有什么关系，在性能上区别大吗. https://www.zhihu.com/question/27872849. 2015/2019-08-16. ↩
-
LiXiang. 比较OpenBLAS，Intel MKL和Eigen的矩阵相乘性能. https://www.leexiang.com/the-performance-of-matrix-multiplication-among-openblas-intel-mkl-and-eigen. 2014-03-14/2019-08-16. ↩
Comments