MATLAB中的GPU编程

多核机器和超线程技术使科学家、工程师和金融分析师能够在各种学科中加速计算密集型应用程序。如今，另一种硬件有望获得更高的计算性能:图形处理单元(GPU)。

gpu最初用于加速图形渲染，现在越来越多地应用于科学计算。与只包含几个核心的传统CPU不同，GPU拥有大量的整数和浮点处理器并行阵列，以及专用的高速内存。一个典型的GPU由数百个这样的小处理器组成(图1)。

然而，GPU大大提高的吞吐量是有代价的。首先，内存访问更有可能成为计算的瓶颈。数据必须在计算前从CPU发送到GPU，然后在计算后从GPU取回。由于GPU通过PCI Express总线连接到主机CPU，因此内存访问速度比传统CPU慢。¹这意味着您的总体计算加速受到算法中发生的数据传输量的限制。其次，用C或Fortran编程gpu需要不同的心智模型和技能集，这可能是困难和耗时的。此外，您必须花时间为特定的GPU微调代码，以优化应用程序的峰值性能。

本文演示了并行计算工具箱中的特性™这使您能够运行MATLAB^®通过对你的代码做一些简单的更改来在GPU上编写代码。我们用谱方法求解一个二阶波动方程来说明这种方法。

为什么要并行化波动方程求解器?

波动方程用于广泛的工程学科，包括地震学、流体动力学、声学和电磁学，用于描述声波、光波和流体波。

使用谱方法求解波动方程的算法是一种很好的并行算法，因为它满足使用GPU加速的两个标准（请参阅“在GPU上执行会加速我的应用程序吗？”）：

它是计算密集型的。该算法执行了许多快速傅里叶变换(fft)和反快速傅里叶变换(ifft)。确切的数字取决于网格的大小(图2)和模拟中包含的时间步长的数量。在不同的矩阵上，每个时间步长需要两个fft和四个ifft，单个计算可能涉及几十万个时间步长。

这是巨大的平行。并行FFT算法被设计为“分而治之”，从而在不同的数据上重复执行类似的任务。此外，该算法需要处理线程之间的大量通信和足够的内存带宽。IFFT可以类似地并行运行。

MATLAB中的GPU计算

在继续介绍波动方程示例之前，让我们快速回顾一下MATLAB如何与GPU一起工作。

FFT、IFFT和线性代数运算是100多个内置MATLAB函数之一，可以通过提供类型为GPUArray的输入参数直接在GPU上执行，GPUArray是并行计算工具箱提供的一种特殊数组类型。这些使能图形处理器的功能是重载换句话说，它们的操作方式取决于传递给它们的参数的数据类型。

例如，下面的代码使用FFT算法在CPU上求一个伪随机数向量的离散傅里叶变换:

A=兰特（2^16,1）；

B=快速傅里叶变换（A）；

为了在GPU上执行相同的操作，我们首先使用gpuArray将数据从MATLAB工作区传输到设备内存的命令。然后我们就可以跑了fft，这是该数据上的重载函数之一:

= gpuArray(兰德(2 ^ 16,1));

B=快速傅里叶变换（A）；

的fft操作是在GPU而不是CPU上执行的，因为它的输入(一个GPUArray)是在GPU上保存的。

结果B存储在GPU上。然而，它在MATLAB工作空间中仍然可见。通过运行B类，我们可以看到它是一个GPUArray。

B类

ans=

parallel.gpu.GPUArray

我们可以继续使用支持GPU的函数在设备上操作B。例如，为了可视化我们的结果情节命令自动在GPUArrays上工作:

情节(B);

要将数据返回到本地MATLAB工作区，可以使用聚集命令例如

C =收集(B);

C现在是MATLAB中的一个双精度对象，可以由任何处理双精度对象的MATLAB函数进行操作。

在这个简单的示例中，执行单个FFT函数所节省的时间通常少于将向量从MATLAB工作区传输到设备内存所花费的时间。这通常是正确的，但取决于硬件和阵列的大小。数据传输开销可能会变得非常大，从而降低应用程序的整体性能，尤其是在CPU和GPU之间反复交换数据以执行相对较少的计算密集型操作的情况下。当数据在GPU上时，对其执行多个操作更有效，仅在需要时才将数据带回CPU²．

请注意，与cpu一样，gpu也有有限的内存。但是，与cpu不同的是，它们不具备与磁盘交换内存的能力。因此，您必须验证您想要保存在GPU上的数据不会超过其内存限制，特别是在处理大型矩阵时。通过运行gpuDevice，您可以查询GPU卡，获取名称、总内存和可用内存等信息。

在MATLAB中实现并加速求解波动方程的算法

为了把上面的例子放到上下文中，让我们在一个真实的问题上实现GPU功能。我们的计算目标是解二阶波动方程

\[\压裂{\部分^ 2 u}{\部分t ^ 2} = \压裂{\部分^ 2 u}{\部分x ^ 2} + \压裂{\部分^ 2 u}{\偏y ^ 2} \]

边界上的条件是u = 0。在空间上采用基于谱法的算法求解，在时间上采用二阶中心差分法求解。

谱方法通常用于求解偏微分方程。使用谱方法，将解近似为连续基函数（如正弦和余弦）的线性组合。在这种情况下，我们采用切比雪夫谱方法，它使用切比雪夫多项式作为基函数。

在每一个时间步，我们用切比雪夫谱方法计算当前解在\(x\)和\(y\)维上的二阶导数。利用这些导数与旧解和当前解结合，采用二阶中心差分法(又称跳蛙法)计算新解。我们选择一个时间步长来保持这种跳变方法的稳定性。

MATLAB算法是计算密集型的，当我们计算解决方案的网格中的元素数量增加时，算法执行所需的时间就会显著增加。当使用2048 x 2048网格在单个CPU上执行时，仅完成50个时间步骤就需要1分钟以上的时间。注意，这一次已经包含了MATLAB中固有的多线程的性能优势。从R2007a开始，MATLAB支持多个函金宝app数的多线程计算。这些函数自动在多个线程上执行，而不需要显式地指定命令来在代码中创建线程。

在考虑如何使用并行计算工具箱加速此计算时，我们将重点关注为每个时间步骤执行计算的代码。图3说明了让算法在GPU上运行所需的更改。注意，计算涉及MATLAB操作，gpu支持的重载函数可通过并行计算工具箱获得。这些操作包括FFT和IFFT、矩阵乘法和各种元素操作。因此，我们不需要以任何方式改变算法来在GPU上执行它。我们简单地将数据传输到GPU使用gpuArray在进入每个时间步计算结果的循环之前。

在GPU上执行计算后，我们将结果从GPU转移到CPU。使能图形处理器的功能引用的每个变量都必须在该图形处理器上创建或转移到该图形处理器上才能使用。

为了将用于光谱微分的权重转换为GPUArray变量，我们使用

W1T = gpuArray (W1T);

某些类型的数组可以直接在GPU上构造，而不需要从MATLAB工作空间转移它们。例如，为了直接在GPU上创建一个零矩阵，我们使用

uxx = parallel.gpu.GPUArray.zeros (N + 1, N + 1);

我们使用聚集功能从GPU带回数据;例如:

vvg=聚集（vv）；

请注意，向GPU进行一次数据传输，然后从GPU进行一次数据传输。每个时间步的所有计算都在GPU上执行。

比较CPU和GPU的执行速度

为了评估使用GPU解决二阶波方程的好处，我们进行了一项基准研究，其中我们测量了在Intel上执行网格大小为64、128、512、1024和2048的50个时间步所花费的时间^®至强^®处理器X5650，然后使用NVIDIA^®特斯拉™C2050 GPU。

对于2048的网格大小，该算法显示计算时间减少了7.5倍，从CPU上的一分钟减少到GPU上的10秒（图4）。对数比例图显示，对于较小的网格大小，CPU实际上更快。然而，随着技术的发展和成熟，GPU解决方案越来越能够处理更小的问题，我们希望这一趋势能够继续下去。金宝搏官方网站

高级图形处理器编程与MATLAB

并行计算工具箱提供了一个直接的方法，通过在GPU上执行它来加快MATLAB代码。您只需更改函数输入的数据类型，以利用许多为GPUArrays重载的MATLAB命令。(支持GPUArray的内建MATLAB函数的完整列表可在金宝app并行计算工具箱文档．）

要在GPU上使用多个简单操作加速算法，可以使用arrayfun，它将函数应用于数组的每个元素。因为arrayfun是一个启用了gpu的功能，您只在单个调用arrayfun，而不是针对每一个单独的操作。

最后，自己编写CUDA代码的有经验的程序员可以使用并行计算工具箱中的CUDAKernel接口将该代码与MATLAB集成。CUDAKernel接口支持更细粒度的控制，以加速那些存在性能瓶颈的代码部分。它创建一个MATLAB对象，该对象提供对已编译成PTX代码(PTX是一个低级并行线程执行指令集)的现有内核的访问。然后调用函数宏指令命令在GPU上评估内核，使用MATLAB数组作为输入和输出。

总结

工程师和科学家成功地使用了GPU技术，该技术原本是为了加速图形渲染，以加速他们的学科特定的计算。只需最少的努力和不需要广泛的gpu知识，您就可以在MATLAB中使用gpu的强大功能。GPUArrays和支持gpu的MATLAB函数帮助您在不使用低级CUDA编程的情况下加快MATLAB操作。如果您已经熟悉gpu编程，那么MATLAB还允许您将现有的CUDA内核集成到MATLAB应用程序中，而不需要任何额外的C编程。

要使用GPU实现加速，您的应用程序必须满足一些标准，其中一个事实是，在CPU和GPU之间发送数据所需的时间必须少于在GPU上运行所获得的性能。如果您的应用程序满足这些标准，那么它是MATLAB提供的GPU功能范围的一个很好的候选者。