在MATLAB中使用gpu

发布的罗兰美国舒尔，2012年2月6日

129次观看(过去30天)|0喜欢|14日的评论

今天我要介绍一位客座博主莎拉等等，扎拉内克他在MATLAB营销团队工作。莎拉之前写关于加快客户的代码以获得可接受的性能。她将讨论如何使用gpu来加速MATLAB应用程序。

内容

概述

在这篇文章中，我们将首先介绍在MATLAB中使用GPU的基础知识，然后使用GPU功能求解二阶波动方程。这篇博客文章的灵感来自于最近的MATLAB文摘文章这是我和我们的一个开发人员吉尔·里斯共同撰写的。从最初的演示开始，MATLAB中可用的GPU函数已经增长。如果您将下面的代码与论文中的代码进行比较，它们略有不同，反映了这些新功能。此外，还做了一些小更改，以便在此博客格式中更容易地解释代码。本文中展示的GPU功能需要并行计算工具箱。

GPU的背景

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

一个好的经验法则是，你的问题可能很适合GPU，如果它是:

大规模并行:计算可以被分解成成百上千个独立的工作单元。当所有核心都处于忙碌状态时，您将看到最佳性能，利用GPU固有的并行特性。看似简单，对含有数十万个元素的数组进行矢量化MATLAB计算通常可以归入这一类。

计算密集型:花费在计算上的时间大大超过了花费在GPU内存之间传输数据的时间。由于GPU通过PCI Express总线连接到主机CPU，因此内存访问比传统CPU慢。这意味着您的整体计算加速受到算法中发生的数据传输量的限制。

不满足这些条件的应用程序实际上可能在GPU上比在CPU上运行得更慢。

了解我们的GPU

有了这些背景知识，我们现在可以开始在MATLAB中使用GPU了。让我们首先查询我们的GPU，看看我们正在处理什么:

gpuDevice

ans = parallel.gpu. cuda设备句柄包:并行。gpu属性:名称:'Tesla C2050 / C2070'索引:1 ComputeCapability: '2.0' SupportsDouble: 1 金宝appDriverVersion: 4 MaxThreadsPerBlock: 1024 MaxShmemPerBlock: 49152 MaxThreadBlockSize: [1024 1024 64] MaxGridSize: [65535 65535] SIMDWidth: 32 TotalMemory: 3.1820e+009 FreeMemory: 2.6005e+009 MultiprocessorCount: 14 ClockRateKHz: 1147000 ComputeMode: 'Default' GPUOverlapsTransfers: 1 KernelExecutionTimeout: 1 CanMapHostMemory: 1 DeviceSupported: 1 DeviceSelected: 1

我们用的是特斯拉C2050。目前，并行计算工具箱支持具有1.3或更高计算能力的NVDIA gpu。金宝app在这里是查看是否支持您的卡的链接。金宝app

一个使用重载函数的简单示例

超过100个操作(例如:fft，传输线，eig)现在可以作为内置的MATLAB函数，通过提供GPUArray类型的输入参数，可以直接在GPU上执行。这些支持gpu的函数是重载的——换句话说，它们的操作取决于传递给它们的参数的数据类型。在这里是所有重载函数的列表。

让我们创建一个GPUArray并执行fft使用GPU。但是，让我们首先在CPU上执行此操作，以便看到代码和性能上的差异

A1 =兰特(3000,3000);抽搐;B1 = fft(A1);Time1 = toc;

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

A2 = gpuArray(A1);抽搐;B2 = fft(A2);Time2 = toc;

第二种情况是在GPU而不是CPU上执行的，因为它的输入(一个GPUArray)保存在GPU上。结果,B2，存储在GPU上。但是，它在MATLAB工作区中仍然可见。通过运行类(B2)，我们可以看到它也是一个GPUArray。

类(B2)

ans = parallel.gpu.GPUArray

为了将数据带回CPU，我们运行收集．

B2 =收集(B2);类(B2)

Ans = double

B2现在在CPU上，并且有一个类double。

加速我们的简单例子

在这个简单的例子中，我们可以计算执行的加速fft根据我们的数据。

speedUp = time1/time2;disp(加速)

3.6122

它看起来像我们的fft在GPU上运行速度快3.5倍。这很好，尤其是自从fft是多线程的核心MATLAB。然而，为了进行真实的比较，我们应该包括将矢量传输到GPU和从GPU传输出去所花费的时间。如果我们这样做，我们会发现加速度大大减小了。让我们看看如果把数据传输的时间也算进去会发生什么。

抽搐;A3 = gpuArray(A1);B3 = fft(A3);B3 =聚(B3);Time3 = toc;speedUp = time1/time3;disp(加速)

0.5676

理解和限制数据传输开销

数据传输开销会变得非常大，以至于会降低应用程序的整体性能，特别是在CPU和GPU之间重复交换数据以执行相对较少的计算密集型操作时。然而，并不是所有的希望都破灭了!通过限制GPU和CPU之间的数据传输，我们仍然可以实现加速。

我们可以直接在GPU上创建数据，而不是在CPU上创建数据并将其传输到GPU上。有几个构造函数可用的以及类构造函数，如meshgrid．让我们看看这对我们的时间有什么影响。为了准确起见，我们还应该重新计算原始串行代码的时间，包括在CPU上生成随机矩阵所需的时间。

抽搐;A4 =兰特(3000,3000);B4 = fft(A4);Time4 = toc;抽搐;A5 = parallel.gpu.GPUArray.rand(3000,3000);B5 = fft(A5);B5 =收集(B5);Time5 = toc;speedUp = time4/time5; disp(speedUp);

1.4100

这是更好的，尽管我们仍然看到从GPU收集数据的影响。在这个简单的演示中，效果被夸大了，因为我们在GPU上运行一个已经快速的单一操作。当数据在GPU上时，对其执行多个操作会更有效，只在需要时才将数据带回CPU。

波动方程的求解

为了把上面的例子放在上下文中，让我们在一个更“真实”的问题上使用相同的GPU功能。为此，我们要解一个二阶波动方程:

求解波动方程的算法结合了空间上的谱法和时间上的二阶中心有限差分法。具体地说，我们应用了切比雪夫谱方法，它使用切比雪夫多项式作为基函数。我们的例子主要是基于Trefethen书中的一个例子:MATLAB中的光谱方法

在GPU上运行算法的代码更改

在加速我们的算法时，我们专注于加速主时间步进while循环中的代码。这部分代码中的操作(例如。fft而且传输线例如，矩阵乘法)都是与GPU一起工作的重载函数。因此，我们不需要以任何方式改变算法来在GPU上执行它。我们可以简单地将数据传输到GPU，完成后再传输回CPU。

类型(“stepSolution.m”）

function [vv,vvold] = stepsolution(vv, vold,ii,N,dt,W1T,W2T,W3T,W4T，…WuxxT1、WuxxT2 WuyyT1 WuyyT2) V = (vv (ii):) vv (ii, N: 1:2)];U = real(fft(V.')).';W1test = (u *W1T).';W2test = (U.*W2T).';W1 = (real(ifft(W1test))).';W2 = (real(ifft(W2test))).';计算x uxx(ii,ii) = W2(:，ii)中的二阶导数。* WuxxT1 - W1(:，ii).*WuxxT2;uxx([1,N+1]，[1,N+1]) = 0;V = [vv(:，ii); vv((N:-1:2),ii)]; U = real(fft(V)); W1 = real(ifft(U.*W3T)); W2 = real(ifft(U.*W4T)); % Calculating 2nd derivative in y uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2; uyy([1,N+1],[1,N+1]) = 0; % Computing new value using 2nd order central finite difference in % time vvnew = 2*vv - vvold + dt*dt*(uxx+uyy); vvold = vv; vv = vvnew; end

将数据传输到GPU并返回的代码更改

的变量dt，x，index1,index2在CPU上进行计算，然后使用gpuArray．

例如:

N = 256;index = 1i*[0:N-1 0 1-N:-1];index1 = gpuArray(index1);

对于其余变量，我们直接在GPU上使用转置运算符的重载版本创建它们(＇)，repmat,meshgrid．

例如:

W1T = repmat(index1, n -1,1);

当我们在GPU上完成所有的计算时，我们使用收集一次将数据从GPU带回来。注意，我们不需要在时间步骤之间将数据传输回CPU。所有这些加在一起看起来是这样的:

类型WaveGPU.m

用谱方法求解二阶波动方程本例求解二阶波动方程:utt = uxx + uyy，边界u = % 0。它在%时间上使用二阶中心有限差分，在空间上使用Chebysehv谱方法(使用FFT)。代码由Trefethen, Lloyd N.在MATLAB中的光谱方法中的一个例子修改，% X和Y中的点X = cos(pi*(0:N)/N);将x发送到GPU x = gpuArray(x);Y = x';计算时间步长dt = 6/N^2;%正在设置网格[xx,yy] = meshgrid(x,y);计算初始值vv = exp(-40*(xx-.4)。^2 + yy.^2));vold = vv;ii = 2:N; index1 = 1i*[0:N-1 0 1-N:-1]; index2 = -[0:N 1-N:-1].^2; % Sending data to the GPU dt = gpuArray(dt); index1 = gpuArray(index1); index2 = gpuArray(index2); % Creating weights used for spectral differentiation W1T = repmat(index1,N-1,1); W2T = repmat(index2,N-1,1); W3T = repmat(index1.',1,N-1); W4T = repmat(index2.',1,N-1); WuxxT1 = repmat((1./(1-x(ii).^2)),N-1,1); WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),N-1,1); WuyyT1 = repmat(1./(1-y(ii).^2),1,N-1); WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,N-1); % Start time-stepping n = 0; while n < Nsteps [vv,vvold] = stepsolution(vv,vvold,ii,N,dt,W1T,W2T,W3T,W4T,... WuxxT1,WuxxT2,WuyyT1,WuyyT2); n = n + 1; end % Gather vvg back from GPU memory when done vvg = gather(vv);