编写可移植的 GPU 代码

本示例演示了如何编写健壮的代码，使其无论在配备 GPU 还是未配备 GPU 的计算机上都能运行。

编写既可在 CPU 上运行，又可在 GPU 上运行的代码，既能让您在不同的计算机上运行该代码并与合作者共享，又能确保任何拥有 GPU 并运行该代码的人都能受益于 GPU 加速。

定义在 CPU 上运行的函数

定义一个名为 matrixOperationsCPU 的函数，该函数执行以下步骤：

计算输入矩阵 A 的逆矩阵
计算逆矩阵的矩阵平方根
创建一个全为 1 的矩阵，其尺寸与输入矩阵 A 相同
将逆矩阵的平方根矩阵与全 1 矩阵相乘

function D = matrixOperationsCPU(A)

% Compute inverse.
B = inv(A);

% Compute matrix square root.
C = sqrtm(B);

% Create matrix of ones with the same size as A.
sz = size(A);
R = ones(sz);

% Multiply matrix square root by matrix of ones.
D = R*C;

end

检查该函数在 CPU 上是否按预期运行。该函数默认在 CPU 上运行。

A = [1 0 2; -1 5 0; 0 3 -9];
C = matrixOperationsCPU(A)

C = 3×3 complex

    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i
    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i
    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i

将编辑函数设置为在 GPU 或 CPU 上运行

有几种选项可以修改该函数，使其既能在 CPU 上运行，也能在 GPU 上运行。最常见的做法是根据输入数据的类型来决定是否使用 GPU。换句话说，如果输入数据是 gpuArray，则在 GPU 上运行该函数；否则，在 CPU 上运行该函数。

有几种函数和语法，可用于创建可在 GPU 或 CPU 上运行的函数。除了 gpuArray 之外，这些功能均无需使用 Parallel Computing Toolbox ™。

gpuArray 和 gather - 使用 gpuArray 函数将数据移动到 GPU 内存中。使用 gather 从 GPU 内存中读取数据。
like 语法--许多数据创建函数（如 rand 和 zeros）采用的 like=p 语法，使得该函数返回的数据与现有数组 p 具有相同的类型和复杂度。
canUseGPU - canUseGPU 函数用于验证是否可用受支持的 GPU，以及 Parallel Computing Toolbox 是否已安装并获得使用许可。
isa - isa 函数用于判断输入数据是否具有指定的数据类型，例如 gpuArray、double 或 single。

如果输入矩阵 A 是 gpuArray，则 matrixOperationsCPU 函数会抛出错误，因为 sqrtm 函数不支持 gpuArray 输入。如果某个函数支持 gpuArray 输入，则会在该函数的参考页面“扩展功能”部分中予以说明。

定义一个新函数 matrixOperations，该函数执行与 matrixOperationsCPU 相同的操作，并包含用于处理 gpuArray 输入的额外代码。该函数：

在使用 sqrtm 函数之前，从 GPU 中读取 B 矩阵。如果 B 不是 gpuArray，那么 gather 函数将不起作用。
创建一个全为 1 的矩阵 R，其类型与输入矩阵 A 相同。如果 A 是 gpuArray，则该函数会直接在 GPU 上创建矩阵 R。

function D = matrixOperations(A)

% Compute inverse.
B = inv(A);

% Compute matrix square root.
B = gather(B);
C = sqrtm(B);

% Create a matrix of ones with the same type and size as A.
sz = size(A);
R = ones(sz,like=A);

% Multiply matrix square root by matrix of ones.
D = R*C;

end

请确认无论是在 CPU 还是 GPU 上运行该函数，其输出结果都是一样的。

C = matrixOperations(A)

C = 3×3 complex

    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i
    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i
    1.1161 - 0.0055i    0.4269 - 0.0555i    0.2283 + 0.2597i

A = gpuArray(A);
CGPU = matrixOperations(A)

CGPU =

   1.1161 - 0.0055i   0.4269 - 0.0555i   0.2283 + 0.2597i
   1.1161 - 0.0055i   0.4269 - 0.0555i   0.2283 + 0.2597i
   1.1161 - 0.0055i   0.4269 - 0.0555i   0.2283 + 0.2597i

请验证当输入为 gpuArray 时，该函数的输出是否为 gpuArray。

isa(CGPU,"gpuArray")

ans = logical
   1

在可移植代码中使用 `gpuArray` 函数

要在可移植代码中使用 gpuArray，请仅在条件语句中调用它，且该语句仅在具备受支持的 GPU 并拥有 Parallel Computing Toolbox 时才会执行。例如，创建一个矩阵，并仅当 canUseGPU 返回 1 (true) 时，才将其转换为 gpuArray。

X = [2 7 9; 2 3 4; 7 1 3];

if canUseGPU
    X = gpuArray(X);
end

定义用于求解二维波方程的函数（选项）

在本节中，您将看到一个更完整、更复杂的代码示例，该代码既可在 CPU 上运行，也可在 GPU 上运行。

waveEquation 函数利用谱方法求解二维波方程 [1]。该函数定义在本示例的末尾，它使用 cast 和 zeros 函数，并采用 like 语法，以确保当输入数据为 gpuArray 时，输出数据为 gpuArray。这使得 waveEquation 函数能够根据输入数据的类型，在 CPU 或 GPU 上运行。

二维波方程可表示为

$\frac{\partial^{2}}{\partial t^{2}} u = \frac{\partial^{2}}{\partial x^{2}} u + \frac{\partial^{2}}{\partial y^{2}} u$ ,

其中 $u = u (x, y, t)$ 是一个标量字段， $x$ 和 $y$ 是空间坐标， $t$ 是时间坐标。该函数假设在二维网格的边界上成立 $u = 0$ 。

指定网格尺寸，并使用切比雪夫点生成网格。

gridSize = 64;

x = cos(pi*(0:gridSize)/gridSize);
y = x';

[xx,yy] = meshgrid(x,y);

定义两个时间步长的边界条件。

vv = exp(-40*((xx-.4).^2 + yy.^2));
vv0 = vv;

在 CPU 上对该方程进行 1000 个时间步长的求解。输入 vv 的类型决定了该函数是在 CPU 还是 GPU 上进行计算。具体来说，如果 vv 是 gpuArray，那么该函数就会在 GPU 上进行计算。

numSteps = 2000;
[vv,xx,yy] = waveEquation(vv,vv0,GridSize=gridSize,NumSteps=numSteps,Plots=true);

比较在 CPU 和 GPU 上运行该函数所需的时间。使用 timeit 函数来计时 CPU 执行，使用 gputimeit 函数来计时 GPU 执行。对于使用 GPU 的函数，gputimeit 比 timeit 更可取，因为它会在记录时间之前检查 GPU 上的所有操作是否已完成，并补偿开销。对于不使用 GPU 的操作，timeit 能提供更高的精度。在 256×256 的网格上运行该函数，共 1000 个时间步长。

gridSize = 256;
numSteps = 1000;
x = cos(pi*(0:gridSize)/gridSize);
y = x';

[xx,yy] = meshgrid(x,y);

vv = exp(-40*((xx-.4).^2 + yy.^2));
vv0 = vv;
tCPU = timeit(@() waveEquation(vv,vv0,GridSize=gridSize,NumSteps=numSteps,Plots=false),3)

tCPU = 
18.7529

vv = gpuArray(vv);
tGPU = gputimeit(@() waveEquation(vv,vv0,GridSize=gridSize,NumSteps=numSteps,Plots=false),3)

tGPU = 
0.8554

disp("Speedup when using GPU: " + round(tCPU/tGPU,1) + "x")

Speedup when using GPU: 21.9x

运行此代码时，性能提升程度将很大程度上取决于您的硬件。

支持函数

`waveEquation` 函数

waveEquation 函数利用谱方法求解二维波方程 [1]。该函数以系统在前两个时间步长时的状态作为输入，并可选地将网格大小、求解的步数以及绘图选项作为输入。

waveEquation 函数采用 cast 和 zeros 函数的语法，并结合 like 的语法，将输入数据指定为原型数组。这使得 waveEquation 函数能够根据输入数据的类型，在 CPU 或 GPU 上以单精度或双精度运行。

function [vv,xx,yy] = waveEquation(vv,vv0,NameValueArgs)

arguments
    % Initial conditions.
    vv  (:,:) {mustBeFloat}
    vv0 (:,:) {mustBeFloat}

    % Name-value arguments.
    NameValueArgs.GridSize (1,1) double {mustBePositive,mustBeInteger} = 64
    NameValueArgs.NumSteps (1,1) double {mustBeNonnegative,mustBeInteger} = 1000
    NameValueArgs.Plots (1,1) logical = true
    NameValueArgs.PlottingFrequency (1,1) double {mustBePositive,mustBeInteger} = 20;
end

% Make GridSize the same type as vv. If vv is a gpuArray, then MATLAB performs subsequent
% operations involving GridSize, including the colon operator, on the GPU.
NameValueArgs.GridSize = cast(NameValueArgs.GridSize,like=vv);

% Set up grid using Chebyshev points.
x = cos(pi*(0:NameValueArgs.GridSize)/NameValueArgs.GridSize);
y = x';
[xx,yy] = meshgrid(x,y);

% Calculate time step.
dt = 6/NameValueArgs.GridSize^2;

% Initialize weights used for spectral differentiation via FFT.
ii = 2:NameValueArgs.GridSize;
index1 = 1i*[0:NameValueArgs.GridSize-1 0 1-NameValueArgs.GridSize:-1];
index2 = -[0:NameValueArgs.GridSize 1-NameValueArgs.GridSize:-1].^2;

W1T = repmat(index1,NameValueArgs.GridSize-1,1);
W2T = repmat(index2,NameValueArgs.GridSize-1,1);
W3T = repmat(index1.',1,NameValueArgs.GridSize-1);
W4T = repmat(index2.',1,NameValueArgs.GridSize-1);

WuxxT1 = repmat((1./(1-x(ii).^2)),NameValueArgs.GridSize-1,1);
WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),NameValueArgs.GridSize-1,1);
WuyyT1 = repmat(1./(1-y(ii).^2),1,NameValueArgs.GridSize-1);
WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,NameValueArgs.GridSize-1);

uxx=zeros(NameValueArgs.GridSize+1,NameValueArgs.GridSize+1,like=vv);
uyy=zeros(NameValueArgs.GridSize+1,NameValueArgs.GridSize+1,like=vv);

for step = 1:NameValueArgs.NumSteps
    V = [vv(ii,:) vv(ii,NameValueArgs.GridSize:-1:2)];
    U = real(fft(V.')).';

    W1test = (U.*W1T).';
    W2test = (U.*W2T).';
    W1 = (real(ifft(W1test))).';
    W2 = (real(ifft(W2test))).';

    % Calculate second derivative in x.
    uxx(ii,ii) = W2(:,ii).* WuxxT1 - W1(:,ii).*WuxxT2;
    uxx([1,NameValueArgs.GridSize+1],[1,NameValueArgs.GridSize+1]) = 0;

    V = [vv(:,ii); vv((NameValueArgs.GridSize:-1:2),ii)];
    U = real(fft(V));

    W1 = real(ifft(U.*W3T));
    W2 = real(ifft(U.*W4T));

    % Calculating second derivative in y
    uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2;
    uyy([1,NameValueArgs.GridSize+1],[1,NameValueArgs.GridSize+1]) = 0;

    % Compute new values using 2nd order central finite difference in time.
    vvnew = 2*vv - vv0 + dt*dt*(uxx+uyy);
    vv0 = vv;
    vv = vvnew;

    % Plot result.
    if NameValueArgs.Plots==true && (step==1 || step==NameValueArgs.NumSteps || mod(step,NameValueArgs.PlottingFrequency)==0)
        updatePlot(xx,yy,vv,step)
    end

end

end

`updatePlot` 函数

waveEquation 函数调用 updatePlot 函数来生成并更新标量场的曲面图。updatePlot 函数接受网格坐标 xx 和 yy、计算得到的解 vv 以及当前步数 step。完成第一步后，该函数会初始化一个曲面图，并设置适当的坐标轴和颜色范围。在后续步骤中，该函数会根据当前值更新图。

function updatePlot(xx,yy,vv,step)

if step == 1
    % Initialize plot.
    figure
    ax = axes;
    surf(ax,xx,yy,vv);
    axis([-1 1 -1 1 -0.1 1]);
    clim([-0.5 0.5]);
    ax.Clipping = "off";
else
    % Update plot.
    ax = gca;
    surface = ax.Children;
    surface.ZData = gather(vv);
end

drawnow
end

参考资料

[1] 特雷费森，洛伊德·N. 《MATLAB》中的谱方法。软件、环境、工具。工业与应用数学学会，2000 年。

另请参阅

gpuArray | gather | canUseGPU | isa