在Matlab上进行CUDA循环

Question

我一直在尝试使用Fortran中的ACC和OpenMP进行并行化。现在我正尝试在Matlab中实现相同的功能。我发现很有趣的是，使用GPU在Matlab中实现循环并行化似乎非常困难。显然，唯一的方法是使用arrayfun函数。但我可能是错的。

从概念上讲，我想知道为什么在Matlab中使用GPU不比在Fortran中更简单。在更实际的层面上，我想知道如何在下面的简单代码中使用GPU。

下面，我分享了三个代码和基准测试：

1. Fortran OpenMP 代码
2. Fortran ACC 代码
3. Matlab parfor 代码
4. Matlab CUDA (?) 这是我不知道如何做的。

Fortran OpenMP 代码：

program rbc

 use omp_lib     ! For timing
 use tools
 implicit none

 real, parameter :: beta = 0.984, eta = 2, alpha = 0.35, delta = 0.01, &
                     rho = 0.95, sigma = 0.005, zmin=-0.0480384, zmax=0.0480384;
integer, parameter :: nz = 4, nk=4800;
real :: zgrid(nz), kgrid(nk), t_tran_z(nz,nz), tran_z(nz,nz);
real :: kmax, kmin, tol, dif, c(nk), r(nk), w(nk);
real, dimension(nk,nz) :: v=0., v0=0., ev=0., c0=0.;
integer :: i, iz, ik, cnt;
logical :: ind(nk);
real(kind=8) :: start, finish   ! For timing
real :: tmpmax, c1  


call omp_set_num_threads(12)


!Grid for productivity z

! [1 x 4] grid of values for z
call linspace(zmin,zmax,nz,zgrid)
zgrid = exp(zgrid)
! [4 x 4] Markov transition matrix of z
tran_z(1,1) = 0.996757
tran_z(1,2) = 0.00324265
tran_z(1,3) = 0
tran_z(1,4) = 0
tran_z(2,1) = 0.000385933
tran_z(2,2) = 0.998441
tran_z(2,3) = 0.00117336
tran_z(2,4) = 0
tran_z(3,1) = 0
tran_z(3,2) = 0.00117336
tran_z(3,3) = 0.998441
tran_z(3,4) = 0.000385933
tran_z(4,1) = 0
tran_z(4,2) = 0
tran_z(4,3) = 0.00324265
tran_z(4,4) = 0.996757

! Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)**(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)**(1/(alpha-1));

! [1 x 4800] grid of possible values of k
call linspace(kmin, kmax, nk, kgrid)


! Compute initial wealth c0(k,z)
do iz=1,nz
  c0(:,iz) = zgrid(iz)*kgrid**alpha + (1-delta)*kgrid;
end do

dif = 10000
tol = 1e-8
cnt = 1

do while(dif>tol)
    !$omp parallel do default(shared) private(ik,iz,i,tmpmax,c1)    
    do ik=1,nk;        
          do iz = 1,nz;
          tmpmax = -huge(0.)

          do i = 1,nk
             c1 = c0(ik,iz) - kgrid(i)
             if(c1<0) exit
             c1 = c1**(1-eta)/(1-eta)+ev(i,iz)
             if(tmpmax<c1) tmpmax = c1
          end do
          v(ik,iz) = tmpmax
       end do

    end do
    !$omp end parallel do
    ev = beta*matmul(v,tran_z)
    dif = maxval(abs(v-v0))
    v0 = v
    if(mod(cnt,1)==0) write(*,*) cnt, ':', dif
        cnt = cnt+1
end do


end program

Fortran ACC:

只需将上述代码中的主循环语法替换为：

do while(dif>tol)
    !$acc kernels
    !$acc loop gang
        do ik=1,nk;        
         !$acc loop gang
          do iz = 1,nz;
          tmpmax = -huge(0.)

          do i = 1,nk
             c1 = c0(ik,iz) - kgrid(i)
             if(c1<0) exit
             c1 = c1**(1-eta)/(1-eta)+ev(i,iz)
             if(tmpmax<c1) tmpmax = c1
          end do
          v(ik,iz) = tmpmax
       end do

    end do

    !$acc end kernels
    ev = beta*matmul(v,tran_z)
    dif = maxval(abs(v-v0))
    v0 = v
    if(mod(cnt,1)==0) write(*,*) cnt, ':', dif
        cnt = cnt+1
end do

Matlab parfor： （我知道下面的代码可以使用向量化语法使其更快，但整个练习的重点是比较循环速度。）

tic;
beta = 0.984; 
eta = 2; 
alpha = 0.35; 
delta = 0.01;
rho = 0.95;
sigma = 0.005;
zmin=-0.0480384;
zmax=0.0480384;
nz = 4;
nk=4800;

v=zeros(nk,nz); 
v0=zeros(nk,nz);
ev=zeros(nk,nz);
c0=zeros(nk,nz);

%Grid for productivity z

%[1 x 4] grid of values for z
zgrid = linspace(zmin,zmax,nz);
zgrid = exp(zgrid);
% [4 x 4] Markov transition matrix of z
tran_z(1,1) = 0.996757;
tran_z(1,2) = 0.00324265;
tran_z(1,3) = 0;
tran_z(1,4) = 0;
tran_z(2,1) = 0.000385933;
tran_z(2,2) = 0.998441;
tran_z(2,3) = 0.00117336;
tran_z(2,4) = 0;
tran_z(3,1) = 0;
tran_z(3,2) = 0.00117336;
tran_z(3,3) = 0.998441;
tran_z(3,4) = 0.000385933;
tran_z(4,1) = 0;
tran_z(4,2) = 0;
tran_z(4,3) = 0.00324265;
tran_z(4,4) = 0.996757;

% Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)^(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)^(1/(alpha-1));

% [1 x 4800] grid of possible values of k
kgrid = linspace(kmin, kmax, nk);

% Compute initial wealth c0(k,z)
for iz=1:nz
  c0(:,iz) = zgrid(iz)*kgrid.^alpha + (1-delta)*kgrid;
end 

dif = 10000;
tol = 1e-8;
cnt = 1;

while dif>tol

    parfor ik=1:nk
          for iz = 1:nz
          tmpmax = -intmax;

          for i = 1:nk
             c1 = c0(ik,iz) - kgrid(i);
             if (c1<0) 
                 continue
             end 
             c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
             if tmpmax<c1 
                 tmpmax = c1;
             end
          end 
          v(ik,iz) = tmpmax;
          end 

    end 
    ev = beta*v*tran_z;
    dif = max(max(abs(v-v0)));
    v0 = v;
    if mod(cnt,1)==0 
        fprintf('%1.5f :  %1.5f \n', [cnt dif])
    end
        cnt = cnt+1;
end 


toc

Matlab CUDA:

这是我不知道如何编码的部分。使用arrayfun是唯一的方法吗？在Fortran中，从OpenMP转换到OpenACC非常简单。在Matlab中，从parfor转换为GPU循环没有简单的方法吗？

代码之间的时间比较：

Fortran OpenMP: 83.1 seconds 
Fortran ACC:    2.4 seconds
Matlab parfor:  1182 seconds

最后一点，我应该说上面的代码解决了一个简单的实际商业周期模型，并基于这个写成。

- phdstudent

---

你的想法是错误的。矢量化语法是MATLAB优化循环的方法 - 无论是在CPU还是GPU上。因此，使用GPU的最简单方法是使用gpuArray()将所有内容放入GPU，然后使用经典的矢量化语法。如果你不能以矢量化方式编写代码，那么arrayfun是一种更繁琐的替代方案。 - Nicky Mattsson

在2018b之前，你不能在MATLAB和CUDA中编写代码。新的、非常专业的工具箱允许在MATLAB中编写内核，但是早期版本只允许使用CUDA运行非常特定的函数。我个人喜欢用CUDA编写代码并将其转换成mex文件。 - Ander Biguri

我正在使用2018b版本。有任何可以查阅的参考资料吗？ - phdstudent

在Matlab中进行性能检查却拒绝使用最基本的优化——向量化是没有意义的。每次迭代都会产生大量开销，而使用向量化操作只会产生一次开销（它本质上也包含了循环）。 - Brice

@phdstudent 是的，官方文档在这里：https://uk.mathworks.com/campaigns/offers/generate-cuda-gpu-code-matlab.html。我认为，除非你需要玩弄主机上分配内存的方式（请参见固定内存），否则最好的方法是编写CUDA代码并将其与Mex混合。 - Ander Biguri

Answer 1

所以，这一部分会给您在这个项目中带来困扰。 MATLAB代表矩阵实验室。向量和矩阵是它的特点。在MATLAB中优化任何东西的主要方法是将其向量化。因此，当使用性能增强工具例如CUDA时，MATLAB假定您会尽可能地向量化输入。鉴于向量化输入在MATLAB编码风格中的重要性，仅使用循环评估其性能并不公平。这就像拒绝使用指针评估C++的性能一样。如果您想在MATLAB中使用CUDA，最主要的方法是向量化输入并使用gpuarray。老实说，我还没有仔细查看您的代码，但它看起来大部分已经向量化了。您可以试试像gpuarray(1:nk)或kgrid=gpuarray(linspace(...)这样简单的方法。

- user10597469

Answer 2

Matlab 编码器

首先，正如Dev-iL已经提到的那样，您可以使用 GPU 编码器。

它（我使用的是 R2019a）只需要对您的代码进行轻微更改：

function cdapted()
beta = 0.984; 
eta = 2; 
alpha = 0.35; 
delta = 0.01;
rho = 0.95;
sigma = 0.005;
zmin=-0.0480384;
zmax=0.0480384;
nz = 4;
nk=4800;

v=zeros(nk,nz); 
v0=zeros(nk,nz);
ev=zeros(nk,nz);
c0=zeros(nk,nz);

%Grid for productivity z

%[1 x 4] grid of values for z
zgrid = linspace(zmin,zmax,nz);
zgrid = exp(zgrid);
% [4 x 4] Markov transition matrix of z
tran_z = zeros([4,4]);
tran_z(1,1) = 0.996757;
tran_z(1,2) = 0.00324265;
tran_z(1,3) = 0;
tran_z(1,4) = 0;
tran_z(2,1) = 0.000385933;
tran_z(2,2) = 0.998441;
tran_z(2,3) = 0.00117336;
tran_z(2,4) = 0;
tran_z(3,1) = 0;
tran_z(3,2) = 0.00117336;
tran_z(3,3) = 0.998441;
tran_z(3,4) = 0.000385933;
tran_z(4,1) = 0;
tran_z(4,2) = 0;
tran_z(4,3) = 0.00324265;
tran_z(4,4) = 0.996757;

% Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)^(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)^(1/(alpha-1));

% [1 x 4800] grid of possible values of k
kgrid = linspace(kmin, kmax, nk);

% Compute initial wealth c0(k,z)
for iz=1:nz
  c0(:,iz) = zgrid(iz)*kgrid.^alpha + (1-delta)*kgrid;
end 

dif = 10000;
tol = 1e-8;
cnt = 1;

while dif>tol

    for ik=1:nk
        for iz = 1:nz
            tmpmax = double(intmin);

            for i = 1:nk
                c1 = c0(ik,iz) - kgrid(i);
                if (c1<0)
                    continue
                end
                c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
                if tmpmax<c1
                    tmpmax = c1;
                end
            end
            v(ik,iz) = tmpmax;
        end

    end
    ev = beta*v*tran_z;
    dif = max(max(abs(v-v0)));
    v0 = v;
    % I've commented out fprintf because double2single cannot handle it 
    % (could be manually uncommented in the converted version if needed)
    % ------------
    % if mod(cnt,1)==0       
    %     fprintf('%1.5f :  %1.5f \n', cnt, dif);
    % end
    cnt = cnt+1;
end
end

构建此内容的脚本如下：

% unload mex files
clear mex


%% Build for gpu, float64

% Produces ".\codegen\mex\cdapted" folder and "cdapted_mex.mexw64"
cfg = coder.gpuConfig('mex');
codegen -config cfg cdapted

% benchmark it (~7.14s on my GTX1080Ti)
timeit(@() cdapted_mex,0)


%% Build for gpu, float32:

% Produces ".\codegen\cdapted\single" folder
scfg = coder.config('single');
codegen -double2single scfg cdapted

% Produces ".\codegen\mex\cdapted_single" folder and "cdapted_single_mex.mexw64"
cfg = coder.gpuConfig('mex');
codegen -config cfg .\codegen\cdapted\single\cdapted_single.m

% benchmark it (~2.09s on my GTX1080Ti)
timeit(@() cdapted_single_mex,0)

因此，如果您的Fortran二进制文件使用float32精度（我怀疑是这样），那么这个Matlab Coder结果与它相当。但这并不意味着两者都非常高效。由Matlab Coder生成的代码仍然远未达到高效的水平。而且它也没有充分利用GPU（即使TDP约为50%）。

---

向量化和 gpuArray

接下来，我同意user10597469和Nicky Mattsson的观点，你的Matlab代码看起来不像正常的“本地”向量化的Matlab代码。

有很多需要调整的地方。（但arrayfun几乎比for好不了多少）。首先，让我们移除for循环：

function vertorized1()
t_tot = tic();
beta = 0.984; 
eta = 2; 
alpha = 0.35; 
delta = 0.01;
rho = 0.95;
sigma = 0.005;
zmin=-0.0480384;
zmax=0.0480384;
nz = 4;
nk=4800;

v=zeros(nk,nz); 
v0=zeros(nk,nz);
ev=zeros(nk,nz);
c0=zeros(nk,nz);

%Grid for productivity z

%[1 x 4] grid of values for z
zgrid = linspace(zmin,zmax,nz);
zgrid = exp(zgrid);
% [4 x 4] Markov transition matrix of z
tran_z = zeros([4,4]);
tran_z(1,1) = 0.996757;
tran_z(1,2) = 0.00324265;
tran_z(1,3) = 0;
tran_z(1,4) = 0;
tran_z(2,1) = 0.000385933;
tran_z(2,2) = 0.998441;
tran_z(2,3) = 0.00117336;
tran_z(2,4) = 0;
tran_z(3,1) = 0;
tran_z(3,2) = 0.00117336;
tran_z(3,3) = 0.998441;
tran_z(3,4) = 0.000385933;
tran_z(4,1) = 0;
tran_z(4,2) = 0;
tran_z(4,3) = 0.00324265;
tran_z(4,4) = 0.996757;

% Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)^(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)^(1/(alpha-1));

% [1 x 4800] grid of possible values of k
kgrid = linspace(kmin, kmax, nk);

% Compute initial wealth c0(k,z)
for iz=1:nz
  c0(:,iz) = zgrid(iz)*kgrid.^alpha + (1-delta)*kgrid;
end 

dif = 10000;
tol = 0.4; 
tol = 1e-8;
cnt = 1;

t_acc=zeros([1,2]);

while dif>tol

    %% orig-noparfor:
    t=tic();
    for ik=1:nk
          for iz = 1:nz
          tmpmax = -intmax;

          for i = 1:nk
             c1 = c0(ik,iz) - kgrid(i);
             if (c1<0) 
                 continue
             end 
             c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
             if tmpmax<c1 
                 tmpmax = c1;
             end
          end 
          v(ik,iz) = tmpmax;
          end 

    end 
    t_acc(1) = t_acc(1) + toc(t);    

    %% better:
    t=tic();          

    kgrid_ = reshape(kgrid,[1 1 numel(kgrid)]);
    c1_ = c0 - kgrid_;    
    c1_x = c1_.^(1-eta)/(1-eta);

    c2 = c1_x + reshape(ev', [1 nz nk]);
    c2(c1_<0) = -Inf;
    v_ = max(c2,[],3);        
    t_acc(2) = t_acc(2) + toc(t);    

    %% compare
    assert(isequal(v_,v));   
    v=v_;

    %% other
    ev = beta*v*tran_z;
    dif = max(max(abs(v-v0)));
    v0 = v;
    if mod(cnt,1)==0
        fprintf('%1.5f :  %1.5f \n', cnt, dif);
    end
    cnt = cnt+1;
end
disp(t_acc);
disp(toc(t_tot));
end

% toc result:
%   tol = 0.4  ->   12 iterations :: t_acc = [  17.7       9.8]
%   tol = 1e-8 -> 1124 iterations :: t_acc = [1758.6     972.0]
% 
% (all 1124 iterations) with commented-out orig :: t_tot = 931.7443

现在，非常明显的是，在while循环内部进行的大多数计算（例如^(1-eta)/(1-eta)）实际上产生了可以预先计算的常量。一旦我们解决了这个问题，结果将比原始的基于parfor的版本（在我的2xE5-2630v3上）快一些：

function vertorized2()
t_tot = tic();
beta = 0.984; 
eta = 2; 
alpha = 0.35; 
delta = 0.01;
rho = 0.95;
sigma = 0.005;
zmin=-0.0480384;
zmax=0.0480384;
nz = 4;
nk=4800;

v=zeros(nk,nz); 
v0=zeros(nk,nz);
ev=zeros(nk,nz);
c0=zeros(nk,nz);

%Grid for productivity z

%[1 x 4] grid of values for z
zgrid = linspace(zmin,zmax,nz);
zgrid = exp(zgrid);
% [4 x 4] Markov transition matrix of z
tran_z = zeros([4,4]);
tran_z(1,1) = 0.996757;
tran_z(1,2) = 0.00324265;
tran_z(1,3) = 0;
tran_z(1,4) = 0;
tran_z(2,1) = 0.000385933;
tran_z(2,2) = 0.998441;
tran_z(2,3) = 0.00117336;
tran_z(2,4) = 0;
tran_z(3,1) = 0;
tran_z(3,2) = 0.00117336;
tran_z(3,3) = 0.998441;
tran_z(3,4) = 0.000385933;
tran_z(4,1) = 0;
tran_z(4,2) = 0;
tran_z(4,3) = 0.00324265;
tran_z(4,4) = 0.996757;

% Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)^(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)^(1/(alpha-1));

% [1 x 4800] grid of possible values of k
kgrid = linspace(kmin, kmax, nk);

% Compute initial wealth c0(k,z)
for iz=1:nz
  c0(:,iz) = zgrid(iz)*kgrid.^alpha + (1-delta)*kgrid;
end 

dif = 10000;
tol = 0.4; 
tol = 1e-8;
cnt = 1;

t_acc=zeros([1,2]);

%% constants:
kgrid_ = reshape(kgrid,[1 1 numel(kgrid)]);
c1_ = c0 - kgrid_;
mask=zeros(size(c1_));
mask(c1_<0)=-Inf;
c1_x = c1_.^(1-eta)/(1-eta);

while dif>tol

    %% orig:
    t=tic();
    parfor ik=1:nk
          for iz = 1:nz
          tmpmax = -intmax;

          for i = 1:nk
             c1 = c0(ik,iz) - kgrid(i);
             if (c1<0) 
                 continue
             end 
             c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
             if tmpmax<c1 
                 tmpmax = c1;
             end
          end 
          v(ik,iz) = tmpmax;
          end 

    end 
    t_acc(1) = t_acc(1) + toc(t);

    %% better:
    t=tic();       
    c2 = c1_x + reshape(ev', [1 nz nk]);
    c2 = c2 + mask;
    v_ = max(c2,[],3);        
    t_acc(2) = t_acc(2) + toc(t);    

    %% compare
    assert(isequal(v_,v));
    v=v_;

    %% other
    ev = beta*v*tran_z;
    dif = max(max(abs(v-v0)));
    v0 = v;
    if mod(cnt,1)==0
        fprintf('%1.5f :  %1.5f \n', cnt, dif);
    end
    cnt = cnt+1;
end
disp(t_acc);
disp(toc(t_tot));
end

% toc result:
%   tol = 0.4  ->   12 iterations :: t_acc = [  2.4    1.7] 
%   tol = 1e-8 -> 1124 iterations :: t_acc = [188.3  115.9]
% 
% (all 1124 iterations) with commented-out orig :: t_tot = 117.6217

这段向量化的代码仍然很低效（例如reshape(ev',...)占用了大约60%的时间，可以通过重新排序维度轻松避免），但它在某种程度上适合于gpuArray()：

function vectorized3g()
t0 = tic();
beta = 0.984; 
eta = 2; 
alpha = 0.35; 
delta = 0.01;
rho = 0.95;
sigma = 0.005;
zmin=-0.0480384;
zmax=0.0480384;
nz = 4;
nk=4800;

v=zeros(nk,nz); 
v0=zeros(nk,nz);
ev=gpuArray(zeros(nk,nz,'single'));
c0=zeros(nk,nz);

%Grid for productivity z

%[1 x 4] grid of values for z
zgrid = linspace(zmin,zmax,nz);
zgrid = exp(zgrid);
% [4 x 4] Markov transition matrix of z
tran_z = zeros([4,4]);
tran_z(1,1) = 0.996757;
tran_z(1,2) = 0.00324265;
tran_z(1,3) = 0;
tran_z(1,4) = 0;
tran_z(2,1) = 0.000385933;
tran_z(2,2) = 0.998441;
tran_z(2,3) = 0.00117336;
tran_z(2,4) = 0;
tran_z(3,1) = 0;
tran_z(3,2) = 0.00117336;
tran_z(3,3) = 0.998441;
tran_z(3,4) = 0.000385933;
tran_z(4,1) = 0;
tran_z(4,2) = 0;
tran_z(4,3) = 0.00324265;
tran_z(4,4) = 0.996757;

% Grid for capital k

kmin = 0.95*(1/(alpha*zgrid(1)))*((1/beta)-1+delta)^(1/(alpha-1));
kmax = 1.05*(1/(alpha*zgrid(nz)))*((1/beta)-1+delta)^(1/(alpha-1));

% [1 x 4800] grid of possible values of k
kgrid = linspace(kmin, kmax, nk);

% Compute initial wealth c0(k,z)
for iz=1:nz
  c0(:,iz) = zgrid(iz)*kgrid.^alpha + (1-delta)*kgrid;
end 

dif = 10000;
tol = 1e-8;
cnt = 1;

t_acc=zeros([1,2]);

%% constants:
kgrid_ = reshape(kgrid,[1 1 numel(kgrid)]);
c1_ = c0 - kgrid_;
mask=gpuArray(zeros(size(c1_),'single'));
mask(c1_<0)=-Inf;
c1_x = c1_.^(1-eta)/(1-eta);

c1_x = gpuArray(single(c1_x));


while dif>tol

    %% orig:
%     t=tic();
%     parfor ik=1:nk
%           for iz = 1:nz
%           tmpmax = -intmax;
% 
%           for i = 1:nk
%              c1 = c0(ik,iz) - kgrid(i);
%              if (c1<0) 
%                  continue
%              end 
%              c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
%              if tmpmax<c1 
%                  tmpmax = c1;
%              end
%           end 
%           v(ik,iz) = tmpmax;
%           end 
% 
%     end 
%     t_acc(1) = t_acc(1) + toc(t);

    %% better:
    t=tic();       
    c2 = c1_x + reshape(ev', [1 nz nk]);
    c2 = c2 + mask;
    v_ = max(c2,[],3);        
    t_acc(2) = t_acc(2) + toc(t);    

    %% compare
    %  assert(isequal(v_,v));        
    v = v_;

    %% other
    ev = beta*v*tran_z;
    dif = max(max(abs(v-v0)));
    v0 = v;
    if mod(cnt,1)==0
        fprintf('%1.5f :  %1.5f \n', cnt, dif);
    end
    cnt = cnt+1;
end
disp(t_acc);
disp(toc(t0));
end

% (all 849 iterations) with commented-out orig :: t_tot = 14.9040

这个 ~15 秒的结果比我们从 Matlab Coder 得到的那些 (~2 秒) 差了 ~7 倍。但是这个选项需要更少的工具箱。在实践中，gpuArray 在你从编写 "本地 Matlab 代码" 开始时最方便。包括交互使用。

最后，如果你用 Matlab Coder 构建这个最终向量化版本 (你需要做一些微不足道的调整)，它不会比第一个更快。它会慢 2-3 倍。

- Igor

---

我接受这个答案，因为它看起来非常棒。我明天稍后会测试一下。谢谢。通过观察代码，我不理解代码如何识别哪个循环可以并行化。或者它们都不能，并且收益在矢量化上？您能发布每个代码运行所需的确切时间吗？ - phdstudent

@phdstudent，好的，我刚刚添加了更多的时间。向量化代码不包含循环（除了while），但底层的矩阵操作库同时使用SIMD和多线程。 - Igor