科学计算 | Matlab 使用 GPU 并行计算

本文转载自： https://sanwen8.cn/p/14bJc10.html

Matlab下直接使用GPU并行计算（预告）<-- 这预告也贴出来太久了，然而我的大论文还是没有写完，但是自己挖的坑一定要填上，我可不是写小说的。

小引言

说它小是因为它只是博士论文的附录一部分，但是其实我还是用了很久才学明白的

中心处理器（CentralProcessing Unit, CPU）是计算机系统的计算和控制核心，在轨道设计中使用计算机程序进行研究和仿真都依赖于CPU的强大计算能力，是研究者经常接触并且熟悉的概念。图形处理器（Graphics Processing Unit, GPU）是计算机设备中用来进行绘图运算，支持显示器设备的芯片。GPU的计算速度和指令复杂度远不及CPU，但是由于其支持高刷新率高分辨率显示设备的需求，GPU具有高并行数、大数据吞吐量的特征，并且在这方面的能力远高于CPU。早在2001年左右，便有人提出了基于GPU的通用计算（General-purpose computing ongraphics processing units, GPGPU）的概念，利用GPU的强大的并行计算能力和相对的低功耗来加速科学计算。基本思路是把需要计算的数据打包成GPU可以处理的图像信息，然后利用处理图像信息的运算来实现科学计算。直到2008年左右，NVIDIA公司推出了CUDA（Compute Unified Device Architecture）并行计算架构，代替了GPGPU的概念，并发布了支持C，C++，Fortran的函数库和编译器。CUDA是NVIDIA发明的一种并行计算平台和编程模型。它通过利用图形处理器 (GPU) 的处理能力，可大幅提升计算性能。此后经过一系列的发展和迭代，现在CUDA已经能够支持符合IEEE标准的双精度运算，并且支持越来越丰富的流程控制，至于是多维数组和显存动态分配等功能。

了解GPU的物理结构和动作方式，在一定程序上有助于提高代码的效率，也有且于开发和调试程序。下面基于本论文使用的GPU型号GTX850M，简要介绍一下GPU的物理结构。如下图所示，该GPU核心的实际运算单位有640个，即每个绿色的小块代表的最基础的计算单元流处理器（Streaming Process，SP），它们在GPU内部被组织成不同的结构，来实现常用的图像处理功能。

上图是计算结构的放大图，一行中的每4个SP和一些共用的存储资源等被组织成一个流多处理器（Stream Multiprocessor，SM）；每8个SM，即每8行，添加一些指令调度功能构成一个warp，是GPU中最小的调度单元。在直接调用CUDA进行C/C++编程时，多个计算thread被组成一个block，然后多个block被组织一个grid，最终GPU按照block来调度所有任务在每个warp上进行计算。可以想象，每个block中thread个数，每个grid中block个数的划分，都会影响到能否充分利用GPU的计算能力，而合理的选择即需要不断的尝试，也需要对并行计算算法的理解和较高水下的编辑能力，是一个极其复杂的过程。

好了，看完这些，实际上并不能帮助你了解GPU大概有什么作用，怎么用，还不如直接去看Matlab帮助文档，上述基本都是废话，然而我相信等你熟悉了用Matlab调用GPU跑程序以后，你会发现这些废话基本还算是精髓。

GPU怎么工具

GPU有显存，还有缓存，也不管它对应上面那个图中哪些块块儿，反正我们也涉及不到这么深入的编程，统称为显存好啦。GPU只能直接计算显存中的数据，就好像CPU只能计算内存中的数据，显存中的数据或者由CPU从内存移入到显存，或者直接由GPU在显存中创建。（很好理解吧，总不能拿个电池在显存上搓一搓就写入了数据吧，你知道我在讲哪个段子对吧，很冷对吧。）

然后，GPU就可以对显存里的数据进行各种计算，但是很显然可以执行哪些计算决定于GPU的能力。比如我的执行 gpuDevice 这个命令后，可以查看显卡的能力 ComputeCapability 是5.0，我只记得1.3以上就可以支持双精度运算，即 double，因此 Matlab 也要求必须是 1.3 以上的 GPU 才可以调用。（一般 gtx 8xx 以上的显卡都是5.0）。

平时不用显卡做数值计算的时候，它其实就是每秒算60几张屏幕大小的图，然后不停地输送给显示器。所以可以想象，显卡的计算能力有多大，然而遗憾的是，喂显示器数据的计算都很简单，简单倒可以把一些算法硬件化来提速，而且单精度就足够了，而且偶尔错一两个像素也无所谓，而且帧数算不过来就跳几帧也看不出来……总之因为它的出身是显卡，这些职业特性导致让它做数值计算会有很多麻烦。当然，Nvidia也在解决，你们肯花钱就好。

好在，用Matlab调用GPU计算，一切都简单了，反正除了修改代码和调研，上述困难我都暂时没有碰到~

在 Matlab 中使用 ode45

这就是本文主要实现的功能，相信有很多人想过用并行计算来加速积分器ode45，但是很遗憾，RKxx这种积分器的数学性质决定它一定是不可以并行的，可以去网上搜，一般中文的答非所问，英文的告诉你不行（有一些其它的积分器有学者在研究并行，基本理念都是用特殊函数基底去近似微分方程的解，然后基函数的系数是有办法可以并行迭代求解的，不展开，外行）。

此处实现的GPU上并行计算ode45是这样，每个GPU的小核计算一个ode45，所以同时可以有1000多个ode45在跑，这个意义上的并行。我的应用场景是搜索参数空间里的格点，每个格点都要积分较长时间，并且在积分中判断每步的状态。

（写得好累，语言的表达总是跟不上思路……）

越写越觉得这个事情不是一句两句能说清楚的，但是实际真的不难……

啰嗦的详细步骤

还是一样的工具箱，在熟悉的parfor下面不远处，有个GPU Computing的条目，建议把它通读一遍，内容很少，不比一篇高水平的文章多。

具体要学**的目录如下

一点都不多，而且只需要学**前三个就好，为了把ode45改写到可以在GPU上运行，我选中的这第三个才要重点理解。

点进去首先看这部分，学**怎么在GPU上建立矩阵

”Transfer Arrays Between Workspace and GPU“基本就是一句语法，

agpu = gpuArray(a)

这样就把正常的一个矩阵a变成了GPU上的矩阵agpu（不用非这么命名，也可以叫axianka，嗯，我是在吐槽）。这是之前讲的CPU把数据从内存移动到显存，可以想象，GPU自动生成数据的功能就是上面那个”Create GPU Arrays Directly“，包含一些随机数函数啥的，而且Matlab一再强调与CPU下的不同，用的时候再看就好了，反正我是用不到。

然后就是是学**怎么用这些gpuArray类型的数据。哦对，上述这些显存上的数据，类型是gpuArray，如下图whos a agpu的结果

想要让CPU用这些gpuArray的时候，用gather可以把它们取回到内存变成double（此处不要问问题，不要多想，我假装啥也不知道，用到时候查查文档或者google一下就可能解决的，相信我~）

继续

然后就是是学**怎么用这些gpuArray类型的数据，学**下图中的Run Built-in Functions on a GPU部分

这部分介绍了Matlab内嵌的，可以支持gpuArray数据类型的数据的函数。就是说，你给这些函数输入之前的agpu，它们是会傻傻地算的，并且计算是发生在GPU上，结果也是在显存中。但是，……算了一会儿下一小段再但是。

再来学**如何在Matlab上运行自定义函数，这里也就是我们将要修改的ode45函数

可以看到，我在这里强调了element-wise函数，在Matlab上运行的自定义函数，必须是每个输入都是1维的，element-wise，就像 .*，./，.^ 什么的一样的element-wise operator。这个element-wise的要求，是对自定义函数内部全部操作而言的，即没有向量操作，没有矩阵操作，没有矩阵乘，没有求逆，没有各种……甚至不能有索引出现，当然，没有前面那些还要索引做什么。我的GPU版本的ode45中矩阵都是折成分量，把乘法一点点写出来的，就跟当年学C时候一样，感觉自己像个小孩子。

哦，值得说一下，在内部不能动态定义数组，比如a = [a,1]这种操作是不可以的，因为CUDA好像不支持这功能，而且，这不就出现数组了嘛，记住原则是element-wise。

然后最令人崩溃的一个表格来了，下图，支持的Matlab code

我当时看到这个的时候，直接略过，这什么嘛，又重复一遍，和前面那个图不是一样样的嘛。

当我在GPU什么都跑不起来的时候，自定义函数不停报错的时候，我又看到这部分的时候，唉，让我静静……再次告诉我，别以为Matlab里有一句是废话，一定仔细看，不然全是坑……

说正经的，这些函数Matlab称之为Supported Matlab Code的意思（我揣测）是说这些函数是所有你可以在GPU上运行的自定义代码里使用的函数，即你写GPU自定义code只能用这些，别的不支持。比如**惯了写成sind(30)的必须现在写成sin(30/180*pi)，明白了？好在pi是支持。仔细看看，其实跟C的基础数学库差不多，多个inf，pi，randi啥的，值得注意的是支持全部的常用流程控制for，while，continue，然而switch，try，catch这些高级货是没有的。

好了，有这些知识，基本上就可以改你手头的代码了，改成符合上述要求的：