Цикл CUDA на Matlab

Итак, этот бит - вот что вас испортит в этом проекте. MATLAB - это аббревиатура от Matrix Laboratory. Векторы и матрицы - это своего рода вещь. Способ номер 1 оптимизировать что-либо в MATLAB - это векторизовать это. По этой причине при использовании инструментов повышения производительности, таких как CUDA, MATLAB предполагает, что вы собираетесь векторизовать свои входные данные, если это возможно. Учитывая первенство векторизации входных данных в стиле кодирования MATLAB, нечестно оценивать его производительность с использованием только циклов. Это было бы похоже на оценку производительности C++ при отказе от использования указателей. Если вы хотите использовать CUDA с MATLAB, основной способ сделать это - векторизовать ваши входные данные и использовать gpuarray. Честно говоря, я не слишком внимательно изучал ваш код, но похоже, что ваши входные данные уже в основном векторизованы. Возможно, вам удастся обойтись чем-то таким простым, как gpuarray(1:nk) или kgrid=gpuarray(linspace(...).

Matlab Coder

Во-первых, как 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, все еще далек от эффективности. И он не полностью использует графический процессор (даже TDP составляет ~ 50%).

Векторизация и gpuArray

Затем я согласен с user10597469 и Ники Мэттссон в том, что ваш код 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 секунд примерно в 7 раз хуже, чем результат (~ 2 секунды), который мы получаем от Matlab Coder. Но для этого варианта требуется меньше наборов инструментов. На практике gpuArray наиболее удобен, когда вы начинаете писать «собственный код Matlab». Включая интерактивное использование.

Наконец, если вы создадите эту окончательную векторизованную версию с помощью Matlab Coder (вам придется внести некоторые тривиальные настройки), она не будет быстрее первой. Было бы в 2–3 раза медленнее.