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

科学计算 | 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这些高级货是没有的。

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

只用了这些基本函数和流程控制

通篇element-wise操作

其它……

（写其它是为了严谨，毕竟我也是新手门外汉二把刀，纯粹是为了提高公众号知名度才来写这**的，才不是因为真心热爱学术喜欢搞技术呢~~话都说这儿了，得帮我推广一下吧，您呢~）

改完以后，用arrayfun这个函数调用，像这样：



arrayfun的第一个参数是在GPU上运行的自定义函数fungpu，这里它内部调用了下面要讲的ode45修改版，后续的参数依次是fungpu的参数，其中只要有一个是gpuArray类型的（这里的t10vector），其它参数全部会被传送到显存，然后fungpu就会在GPU上运行。

下面我就开始介绍修改ode45的过程，

（其实想想你跟我修改ode45并没有多大用，重点还是上面讲的那些多看看多琢磨，平时科研碰到用得上的问题留个心，有功夫了就尝试一下，有问题多思考或者去google或者来交流讨论，这样才有效果）

考虑到这么深切的原因而不是因为我懒，我就精简地讲一下。

先看下ode45的源码，我知道其实很多人都已经读过了。

没读过也不要紧，相信我，如果不是有需求被逼，你情愿自己写一个rk45或者rk78也不愿意读这功能强大的代码（不得不说官方的代码写得很好，好用又好读，还好改）





上面红杠的部分，基本是要删掉的，因为在GPU上不支持，这些都是一些高级功能，比如Event啊啥的，还有一些通用性的功能。

我们需要的功能主要在红框中，主要是ode45的rk45算法部分。

（注意上面的代码我折叠了很多，折叠了很多，折叠了很多，而且不是严格按照我划的改，随手一划，你懂得，随手）

修改这部分代码

比如我是积分三维空间中的轨道，所以输入的变量是6维，时间是1维，还有些其它质量比例什么的参数。

改后的结果如下：可以看到此时精度也作为独标量直接输入，而不能再用Name-Value pair的形式输入，其它的功能都被精简掉了，因为GPU不支持。



可以看到原先的矩阵操作都被写成了分量形式，其实也不难，对吧。

下面看到的一些中文注释，很多是在读ode45代码时根据当时理解添加的。注释总不嫌多，对吧。



看，RK45每步求微分的过程，也要用分量形式写出来了，其中的OdeRtbpElementWise函数应该能猜到是什么吧，就是要积分的动力学方程，也被写成了分量形式的输入输出，1个时间，6个状态。



误差分析也要一点一点拆，就当复**了一遍数值分析，对吧。你看这大段大段的英文原版注释，我都没有动过的，说明修改的工作量其实也不大，对吧。



这就结束了，有没有觉得很神奇？这货就可以在GPU上积分啦！

等等，并没有，还没有讲输出。因为GPU不支持自定义函数的索引功能，不支持动态数组，我们不能像正常ode45一样输出轨道的时刻序列，只能输出某一个状态，这里我们选择末状态[t,y1,y2,y3,y4,y5,y6]。

这样便存在一个问题，如何验证整个过程中轨道是正确的？很简单。

用ode45从t0到t1积分一遍，取ii=1:1000个点画图。

用ode45gpu积分一遍，分别从t0到t0+ii/1000积分ii=1:1000条轨道，注意这1000条轨道是同时被计算的，所以和从t0到t1直接积分的耗时只差一点点一点点。

然后gather回来重叠画图比较


这个图就是这么画的，所以你看这么多小点点~~怎么样，不难吧？

还有一个调试的问题没有讲到，是程序就免不了调试。其实也不用讲，都在Matlab下了，还要怎么讲，直接在CPU下串行调试就好啦，保证element-wise的原则，正常调试喽。

应该就讲完了，再讲就是各种细节了，我觉得没必要。

最后上个图说明一下有效性：



上图是在星历模型下搜索得到的某种结果图，下图是在某种近似模型下搜索得到的同种结果图，关键点在于分布规律是一样的！说明GPU的积分精度和直接在CPU上的RK78积分精度是可比拟的，是可用的。

上图大约14小时，分12x16线程计算，2.3GHz的CPU，用的UPC的服务器，还是FORTRAN77的计算效率；

下图大约13小时，显卡GTX850M超频后1.2GHz，Matlab自动划分各种Warp，Block啥的。

可见，加速能力还是相当强悍的，才是个不到1000块钱的笔记本显卡。

当然这里有一个问题，GPU如何使用星历模型，窃以为这个问题并不存在一个简单的可以仅依靠Matlab解决的途径，所以，管它呢……

写得着实太粗糙，如果有问题或者有错误，还请右下角随意“留言”，我会认真回复的。

如果有想以我这个为基础研究一下的，我可以考虑分享我的代码给你，暂时还没有想好形式，而且应该要在6月份以后了。

Matlab的SVN版本控制学**（R2014b）

如此这般优雅地TypeMath

原创插画 | 你恰好也在

工具心得 | 未经严谨验证的 Word 崩溃救命技巧

欢迎留言，欢迎讨论。

欢迎推荐给其它人。

点击“阅读原文”查看更多。