利用硬件加速大规模科学计算 | 从MATLAB到C/C++代码
作者 周拥华、MathWorks
无论是传统的系统仿真还是如今火热的人工智能，都涉及到大量的科学计算，高手们也许对利用多处理器、集群和GPU的编程技术驾轻就熟，但对初学者而言怎么样利用硬件来加速大规模科学计算无疑是个门槛较高的问题。作者试图用这篇简单的文章，帮助初学者了解，从MATLAB到C/C++代码，可供使用的一些简单快捷的选项，算是入门吧。
1. 用CPU做矩阵基本运算
假如我有两个10000x10000的矩阵（没错，1亿个元素），A和B，要对它俩进行相乘，得到C，用MATLAB表述是这样的：
A = rand(10000, 10000); % 随机产生10000x10000的矩阵
B = rand(10000, 10000);
C = A*B;
如果简单直白地用C来实现矩阵乘法，它大概是这样子的：
static double a[100000000]; static double b[100000000]; static double c[100000000]; double d; int i; int k; int xpageoffset; for (i = 0; i < 10000; i++) { for (xpageoffset = 0; xpageoffset < 10000; xpageoffset++) { d = 0.0; for (k = 0; k < 10000; k++) { d += a[i + 10000 * k] * b[k + 10000 * xpageoffset]; } c[i + 10000 * xpageoffset] = d; } }
事实上，上述代码正是我从MATLAB自动生成的C代码里删节出来的。
将一个MATLAB里编写的函数或脚本文件生成C代码很简单，你可以通过APP菜单里的MATLAB Coder按提示一步一步来做，也可以通过命令行来实现，譬如下面几行指令可以将一个名为largeMatrixTest.m的脚本文件转换成C代码，并编译为exe（借助MinGW或Visual C++）：
cfg = coder.config('exe'); cfg.GenerateExampleMain = 'GenerateCodeAndCompile'; codegen largeMatrixTest -config cfg -report
如果你在同一台计算机上（譬如我的Lenovo T480s），分别执行第一部分的MATLAB脚本和编译执行第二部分的C代码，效率会差多少呢？我得到的分别是20.5秒和……超过30分钟吧，真的没耐心等它执行完。那么问题来了，为啥MATLAB会比C还快？
如果你关注一下两者运行时CPU的占用率，你会发现，前者占用了你CPU的大多数核，CPU总占用率达到了70%以上，而后者在我4核8线程的机器上只占用了12%的CPU。MATLAB是怎么做到的？答案是BLAS（Basic Linear Algebra Subprograms），它是一个为底层向量与矩阵运算针对具体处理器高度优化实现的库，通过BLAS，大规模矩阵运算便能利用计算机的多核来加速。BLAS有C接口，我们当然可以手写C代码实现同样的加速，但显然没MATLAB那么容易.
有没有可能让MATLAB生成的C代码也能用到BLAS呢？答案是肯定的，参考这里就好: MATLAB帮助文档Speed Up Matrix Operations in Generated Standalone Code by Using BLAS Calls
参照上面的文档，从MATLAB生成的C代码文件中我们可以看到，矩阵乘法运算C=A*B对应的代码是这样的：
static double a[100000000]; static double b[100000000]; static double c[100000000]; cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, (MKL_INT)10000, (MKL_INT)10000, (MKL_INT)10000, 1.0, &a[0], (MKL_INT)10000, &b[0], (MKL_INT)10000, 0.0, &c[0], (MKL_INT)10000);
通过代码生成的报告我们还可以了解m脚本和生成的C代码之间的对应关系（如本文开头的图片所示）。
经过BLAS加持的C代码，跑起来效率就跟MATLAB脚本一样快了。
2. 用CPU做矩阵求解
BLAS就够用了吗？不够。BLAS实现了一些底层的计算，但如果要对矩阵进行求解（如奇异值分解或求特征值），我们还得借助LAPACK来加速。
MATLAB 中的 LAPACK
LAPACK（线性代数包）是一个例程库，它为数值线性代数和矩阵计算提供快速、稳健的算法。MATLAB 中的线性代数函数和矩阵运算均基于 LAPACK 构建，并且继续受益于其例程的性能和精度。
简史
MATLAB 诞生于 20 世纪 70 年代后期，是一款基于 LINPACK 和 EISPACK 构建的交互式计算器，而 LINPACK 和 EISPACK 在当时是进行矩阵计算的最先进的 Fortran 子例程库。多年来，MATLAB 使用了 LINPACK 和 EISPACK 的十几个 Fortran 子例程的 C 语言版本。
2000 年，MATLAB 改用 LAPACK，这是 LINPACK 和 EISPACK 的现代替代品。它是一个用于数值线性代数的大型、多作者 Fortran 库。LAPACK 最初是为在超级计算机上使用而设计的，因为它能够一次计算矩阵的多个列。LAPACK 例程的速度与基本线性代数子例程 (BLAS) 的速度密切相关。BLAS 版本通常特定于硬件并经过高度优化。
详见：https://www.mathworks.com/company/newsletters/articles/matlab-incorporates-lapack.html
参考这个链接，可将MATLAB脚本中有关矩阵求解的函数（如svd、eig等）生成LAPACK调用。
3. 为for循环加速
看了前面两节的内容，我们已经知道，在MATLAB中，能直接对向量或矩阵做基本运算或调用某个函数去求解就不要自己写for循环来做了，因为MATLAB的向量/矩阵计算和求解基本都是BLAS和LAPACK加持过的。
但我们在实际计算中还是不可避免要有自己写的for循环，有些是顶层的循环，还有些可能是底层的没有可以直接拿来用的函数就自己写了。那么问题来了，自己m脚本里的for循环，无论是在MATLAB里还是生成C代码，跑起来都觉得慢怎么办？
我们知道，在MATLAB里把for换成parfor，就可以通过并行计算工具箱利用单机多处理器或计算机集群硬件来加速，例如：
function a = test_parfor %#codegen a=ones(10,256); r=rand(10,256); parfor (i=1:10, 4) % 指定4线程并发 a(i,:)=real(fft(r(i,:))); end % codegen -config:lib test_parfor
那如果对上述脚本生成C代码，还能利用单机多核来加速吗？可以！
我们也可以从生成的代码中看到（截取部分如下所示），parfor被实现成了多线程（用到了OpenMP，详情可查阅parfor的帮助文档）。我们还可以通过parfor的参数来指定线程的数量，如parfor (i=1:10, 4)，即指定为4线程。
#pragma omp parallel for \ num_threads(4 > omp_get_max_threads() ? omp_get_max_threads() : 4) \ private(b_i,yCol,b_r) for (i = 0; i < 10; i++) { /* 指定4线程并发 */ for (b_i = 0; b_i < 256; b_i++) { b_r[b_i] = r[i + 10 * b_i]; } c_FFTImplementationCallback_doH(b_r, 0, yCol); for (b_i = 0; b_i < 256; b_i++) { a[i + 10 * b_i] = yCol[b_i].re; } }
当然，如果你机器的CPU核数有限或者内存有限，尤其在BLAS和LAPACK或FFTW的并行化已经占用了大多数资源的情况下，期望再通过顶层的parfor提高整体运算效率可能就比较难了，此时你得先提高机器配置，让BLAS、LAPACK或FFTW能尽情施展后，再通过parfor激活for循环的并行化。又或者你的for循环里没有可通过BLAS、LAPACK或FFTW实现的并行化或者计算负载主要在非并行化的部分，就赶紧参考这个链接把parfor用起来吧。
4. 用GPU加速
大家可能听说过深度学习往往要用GPU来加速，这是因为深度学习的模型训练和实时推断都有超大的计算量。先撇开深度学习，在我们经典的科学计算中，怎么用GPU来加速呢？
首先，MATLAB中利用GPU的方法很简单，我们来看下面的例子，C = A*B然后对C矩阵做奇异值分解，你需要做的只不过是把矩阵A和B放进GPU内存罢了，MATLAB会自动将后续的矩阵相乘和奇异值分解分配到GPU上完成。看代码：
function s = largeMatrixTest() coder.gpu.kernelfun; a = rand(5000, 5000); b = rand(5000, 5000); %c = a * b; gpuA = gpuArray(a); gpuB = gpuArray(b); c = gpuA * gpuB; s = svd(c); end % 执行下面的指令，可以统计运算所耗时间（与CPU上不同，用GPU做计算要用wait）： dev=gpuDevice();tim=tic();largeMatrixTest;wait(dev);gpuTime=toc(tim);
尽管多出了两个内存搬运的操作，利用NVIDIA GeForce GTX 1080运行上面的脚本，结果比我T480s的CPU上跑相同的操作要38秒还是快了很多，只需要11秒了。那么大家应该知道，这归功于CUDA。
上面的脚本能直接生成CUDA代码吗（C++）？然而并不能，因为gpuArray并不支持代码生成。那怎么才能生成调用CUDA库的C++代码呢？答案是用最初的针对CPU的脚本，在生成时代码时指定用GPU Coder就行了。来看m代码：
%% largeMatrixTest.m function s = largeMatrixTest() coder.gpu.kernelfun; a = rand(5000, 5000); b = rand(5000, 5000); c = a * b; %gpuA = gpuArray(a); %gpuB = gpuArray(b); %c = gpuA * gpuB; s = svd(c); end %% generate standalone exe by using GPU Coder (gpuCodeGenTest.m) cfg = coder.gpuConfig('exe'); cfg.GenerateExampleMain = 'GenerateCodeAndCompile'; codegen largeMatrixTest -config cfg -report
是不是太简单了，只是将CPU代码生成时的coder.Config换成了coder.gpuConfig，没错，就这么简单！但生成的代码就复杂多了，限于篇幅，这里就只截取C=A*B对应的函数和svd奇异值分解对应的函数，分别如下：
cublasDgemm(getCublasGlobalHandle(), CUBLAS_OP_N, CUBLAS_OP_N, 5000, 5000,
5000, (double *)gpu_alpha1, (double *)&(*gpu_a)[0], 5000, (double *)
&(*gpu_b)[0], 5000, (double *)gpu_beta1, (double *)&(*gpu_c)[0],
5000);
cusolverDnDgesvd(getCuSolverGlobalHandle(), 'N', 'N', 5000, 5000, (double *)
&(*gpu_c)[0], 5000, &(*gpu_s)[0], NULL, 1, NULL, 1, (double
*)getCuSolverWorkspaceBuff(), *getCuSolverWorkspaceReq(), &(*gpu_superb)[0],
gpu_info_t);
显然，这里分别用到了CUDA的cuBLAS和cuSOLVER，另外，CUDA也有cuFFT。
如果你看的仔细，你可能还注意到了在largeMatrixTest.m这个脚本中，有一行特别的代码，coder.gpu.kernelfun，这是一行不影响执行但会影响代码生成的脚本，它告诉GPU Coder，在为这个函数生成C++代码时，将计算任务尽可能映射到GPU上去，MATLAB会对脚本中的for循环进行分析，在允许的情况下为其生成为可并行化执行的函数。除此之位，我们还可以显式地通过coder.gpu.kernel将紧随其后的for循环加以并行化。
【备注】截至R2020b，MATLAB共提供705个可通过gpuArray利用NVIDIA GPU进行加速计算的函数；有326个函数可以生成利用GPU加速的C++代码（即CUDA代码）。除此之外，MATLAB内嵌并行化支持的函数也有62个，覆盖了深度学习训练、图像的非重复块处理、非线性优化、统计与机器学习、特征提取与工业统计等方方面面。众所周知，在计算性能方面，近20年来MATLAB的几乎每个版本更新都有所增强，如果你关心计算性能，用最新的版本是最好的选择。
5. 能支持ARM么？
上面的例子，无论是Intel MKL还是CUDA，默认都是x86/x86_64平台的，操作系统可能是Windows，也可能是Linux，我不确切它们在mac OS会怎样，也许你可以自己尝试一下。
如果你在用基于ARM的系统，可以肯定的一点是，MATLAB目前并没有能运行在ARM上的版本，也许以后会有，但如果你想把MATLAB上开发的算法生成C/C++代码，跑到ARM上可以么？答案是肯定的。那如果要利用ARM上的特定核来加速行不行呢？
如果你的ARM上有GPU并支持CUDA，那么可以肯定，上述第4节的内容依然适用，事实上，MATLAB现在可以非常好地支持针对NVIDIA Jetson和Drive两个系列的产品的代码生成（参见GPU Coder相关文档）。
那如果没有GPU，仅用ARM处理器能生成并行化加速的C代码么？答案应该是ARM+Linux架构上的OpenBLAS、CLAPACK和FFT，以及OpenMP基础库。此外，ARM也提供了一个Arm Performance Libraries (commercial variant) ，大家可以去尝试一下。
【注】作者没有对ARM上的BLAS、LAPACK和FFT做任何测试，理论上可行，实际使用中如遇到问题，建议询问MathWorks中国办公室获得技术支持。
6. 其它
除了以上提到的内容，如今最热且重度依赖硬件加速的深度学习应用并没在本文中讨论，事实上MATLAB从R2017b就已经开始支持针对深度学习推断生成C/C++代码，并可利用硬件来加速深度学习的推断，包括NVIDIA的桌面与服务器GPU及嵌入式GPU（通过CUDA实现）、ARM Mali GPU与ARM Neon核（通过Arm Compute Library实现），或者利用x86_64处理器的SIMD（SSE/AVX，通过Intel MKL-DNN实现）。在最新的R2020b版本中，Deep Learning HDL Toolbox还可以将训练好的深度学习模型生成为硬件描述语言，从而把深度学习部署到FPGA上。详情可参考MATLAB帮助文档或者咨询MathWorks中国办公室。
附录：利用BLAS和LAPACK的完整示例脚本
要从MATLAB生成C代码时支持BLAS和LAPACK，你需要参照前面有关链接安装INTEL MKL（含BLAS和LAPACK），并从coder.BLASCallback和coder.LAPACKCallback分别派生自己的类，用来指定编译所需的头文件、库文件，它们的路径，以及所需的宏定义，最后在代码生成时指定这两个类，样例脚本如下：
%% example function for code generation (largeMatrixTest.m) function largeMatrixTest() a = rand(5000, 5000); tic; b = a * a; c = sum(a); s = svd(a); e = eig(a); [maxValue, maxPos] = max(a); tCpu = toc; fprintf(' Time cost: %f\n', tCpu); end %% define class for BLASCallback (useMyBLAS.m) classdef useMyBLAS < coder.BLASCallback methods (Static) function updateBuildInfo(buildInfo, ~) libPath = 'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows\mkl\lib\intel64'; libPriority = ''; libPreCompiled = true; libLinkOnly = true; libs = {'mkl_intel_ilp64.lib' 'mkl_intel_thread.lib' 'mkl_core.lib'}; buildInfo.addLinkObjects(libs, libPath, libPriority, libPreCompiled, ... libLinkOnly); buildInfo.addLinkObjects('libiomp5md.lib',fullfile(matlabroot,'bin', ... 'win64'), libPriority, libPreCompiled, libLinkOnly); buildInfo.addIncludePaths('C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2020.1.216\windows\mkl\include'); buildInfo.addDefines('-DMKL_ILP64'); end function headerName = getHeaderFilename() headerName = 'mkl_cblas.h'; end function intTypeName = getBLASIntTypeName() intTypeName = 'MKL_INT'; end function doubleComplexTypeName = getBLASDoubleComplexTypeName() doubleComplexTypeName = 'my_double_complex_type'; end function singleComplexTypeName = getBLASSingleComplexTypeName() singleComplexTypeName = 'my_single_complex_type'; end function p = useEnumNameRatherThanTypedef() p = true; end end end %% define class for LAPACKCallback (useMyLAPACK.m) classdef useMyLAPACK < coder.LAPACKCallback methods (Static) function hn = getHeaderFilename() hn = 'mkl_lapacke.h'; end function updateBuildInfo(buildInfo, buildctx) buildInfo.addIncludePaths(fullfile(pwd,'include')); libName = 'mkl_lapack95_ilp64'; libPath = 'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows\mkl\lib\intel64'; [~,linkLibExt] = buildctx.getStdLibInfo(); buildInfo.addLinkObjects([libName linkLibExt], libPath, ... '', true, true); buildInfo.addIncludePaths('C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2020.1.216\windows\mkl\include'); buildInfo.addDefines('HAVE_LAPACK_CONFIG_H'); buildInfo.addDefines('LAPACK_COMPLEX_STRUCTURE'); buildInfo.addDefines('LAPACK_ILP64'); end end end %% generate standalone exe for above MATLAB function (genCodeTest.m) cfg = coder.config('exe'); cfg.CustomBLASCallback = 'useMyBLAS'; cfg.CustomLAPACKCallback = 'useMyLAPACK'; cfg.GenerateExampleMain = 'GenerateCodeAndCompile'; codegen largeMatrixTest -config cfg -report
2022 年发布
