个人中心
gpt4 book ai didi

matlab - MATLAB 上的 CUDA 循环

转载 作者：太空宇宙 更新时间：2023-11-03 19:25:42 24 4
gpt4 key购买 nike

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

在概念层面上，我想知道为什么 GPU 在 matlab 中的使用并不比在 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 参数:(我知道使用矢量化语法可以使下面的代码更快，但练习的重点是比较循环速度)。

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

最后，我应该说上面的代码解决了一个简单的真实商业周期模型，并且是基于 this 编写的.

最佳答案

Matlab 编码器

首先，作为Dev-iL already mentioned ，您可以使用 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 代码看起来不像普通的“ native ”矢量化 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 倍。但是这个选项需要更少的工具箱。实际上，当您开始编写“原生 Matlab 代码”时，gpuArray 最方便。包括交互式使用。

最后，如果您使用 Matlab Coder 构建这个最终的矢量化版本(您将不得不做一些微不足道的调整)，它不会比第一个更快。它会慢 2 到 3 倍。

关于matlab - MATLAB 上的 CUDA 循环，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/53374053/

24 4 0
文章推荐： c# - 如何在前瞻中匹配字符串的结尾？
文章推荐： html - CSS 样式表无法使用 Ruby on Rails:Sprockets::Rails::Helper::AssetNotFound 错误
文章推荐： html - 在 CSS 中，如何使一个元素相对于包含元素居中，而不考虑其他同级元素？
文章推荐： c# - FileStream Write() 方法仅写入未知符号
太空宇宙
个人简介

我是一名优秀的程序员,十分优秀！

滴滴打车优惠券免费领取
滴滴打车优惠券
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com