时间戳

首页 > 解决方案 > matlab上的CUDA循环

matlab - matlab上的CUDA循环

问题描述

我一直在使用 Fortran 中的 ACC 和 OpenMP 进行并行化。我现在正在尝试在 matlab 中做同样的事情。我觉得非常有趣的是，在 matlab 中使用 GPU 并行化循环似乎非常困难。显然，唯一的方法是使用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 加速器：

只需将上面代码中的 mainloop 语法替换为：

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

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

标签： matlabopenmpopenacc

解决方案

所以，这一点会让你在这个项目上搞砸。MATLAB 代表矩阵实验室。向量和矩阵就是它的东西。在 MATLAB 中优化任何东西的第 1 种方法是对其进行矢量化。因此，在使用 CUDA 等性能增强工具时，MATLAB 假定您将尽可能对输入进行矢量化处理。鉴于在 MATLAB 编码风格中向量化输入的重要性，仅使用循环来评估其性能并不是一个公平的比较。这就像在拒绝使用指针的同时评估 C++ 的性能。如果你想在 MATLAB 中使用 CUDA，主要的方法是将你的输入向量化并使用 gpuarray。老实说，我并没有仔细研究您的代码，但看起来您的输入已经大部分是矢量化的。gpuarray(1:nk)或kgrid=gpuarray(linspace(...)。

推荐阅读