MATLAB图形处理器编程

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

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

然而，GPU极大地提高了吞吐量是有代价的。首先，内存访问更有可能成为计算的瓶颈。在计算之前，数据必须从CPU发送到GPU，然后再从GPU中检索。由于GPU通过PCI Express总线连接到主机CPU，因此内存访问比传统CPU慢。¹这意味着您的整体计算加速受到算法中发生的数据传输量的限制。其次，用C或Fortran为图形处理器编程需要不同的思维模式和技能集，获得这些技能很难且耗时。此外，您必须花时间针对特定的GPU对代码进行微调，以优化应用程序以获得最佳性能。

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

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

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

使用频谱方法求解波动方程的算法是并行化的一个很好的候选者，因为它满足使用GPU加速的两个标准(参见“GPU上的执行会加速我的应用程序吗?”):

这需要大量的计算。该算法执行多次快速傅里叶变换(FFTs)和快速傅里叶反变换(IFFTs)。确切的数字取决于网格的大小(图2)和模拟中包含的时间步骤的数量。每个时间步需要在不同的矩阵上进行两个fft和四个ifft，一次计算可能涉及数十万个时间步。

它是高度平行的。并行FFT算法旨在“分而治之”，以便在不同的数据上重复执行类似的任务。此外，该算法需要处理线程之间的大量通信和大量的内存带宽。IFFT也可以并行运行。

MATLAB中的GPU计算

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

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

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

A = rand(2^16,1);

B = fft(A);

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

A = gpuArray(rand(2^16,1));

B = fft(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²．

注意，gpu和cpu一样，内存是有限的。但是，与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个时间步骤就需要超过一分钟。请注意，这一次已经包括了MATLAB中固有多线程的性能优势。从R2007a开始，MATLAB支持多个函万博1manbetx数的多线程计算。这些函数自动在多个线程上执行，而不需要显式地指定命令在代码中创建线程。

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

在GPU上执行计算后，我们将结果从GPU传输到CPU。GPU功能所引用的变量必须在GPU上创建或传输到GPU上才能使用。

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

W1T = gpuArray(W1T);

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

uxx = parallel.gpu. gpuarray .zero (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代码来加速MATLAB代码的直接方法。您只需更改函数输入的数据类型，以利用为gpu阵列重载的许多MATLAB命令。(支持GPUArray的内置MATLAB函数的完整列表可在万博1manbetx并行计算工具箱文档)。

要在GPU上使用多个简单操作加速算法，可以使用arrayfun，它对数组的每个元素应用一个函数。因为arrayfun是一个gpu启用的功能，您只在单个调用上招致内存传输开销arrayfun，而不是在每个单独的操作上。

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

总结

工程师和科学家正在成功地使用GPU技术，该技术最初用于加速图形渲染，以加速他们特定学科的计算。只需很少的努力，并且没有广泛的gpu知识，您现在就可以在MATLAB中使用gpu的强大功能。gpu阵列和gpu支持的MATLAB函数帮助您加快MATLAB操作，而无需低级CUDA编程。如果您已经熟悉gpu编程，MATLAB还可以让您将现有的CUDA内核集成到MATLAB应用程序中，而不需要任何额外的C编程。

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