title: CUDA loop on matlab
body: <p>I have been playing around with parallelization both using ACC and OpenMP in Fortran. I am now trying to do the same in matlab. I find it very interesting that it seems to be very hard to paralelize a loop using GPUs in matlab. Apparently the only way to do it is to by using <code>arrayfun</code> function. But I might be wrong. </p>
  
  <p>At a conceptual level, I am wondering why is the GPU usage in matlab not more straightforward than in fortran. At a more practical level, I am wondering how to use GPUs on the simple code below. </p>
  
  <p>Below, I am sharing three codes and benchmarks:</p>
  
  <ol>
  <li>Fortran OpenMP code</li>
  <li>Fortran ACC code</li>
  <li>Matlab parfor code</li>
  <li>Matlab CUDA (?) this is the one I don't know how to do.</li>
  </ol>
  
  <p><strong>Fortran OpenMP:</strong></p>
  
  <pre><code>program rbc
  
   use omp_lib     ! For timing
   use tools
   implicit none
  
   real, parameter :: beta = 0.984, eta = 2, alpha = 0.35, delta = 0.01, &amp;
                       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&gt;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&lt;0) exit
               c1 = c1**(1-eta)/(1-eta)+ev(i,iz)
               if(tmpmax&lt;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
  </code></pre>
  
  <p><strong>Fortran ACC:</strong></p>
  
  <p>Just replace the mainloop syntax on the above code with:</p>
  
  <pre><code>do while(dif&gt;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&lt;0) exit
               c1 = c1**(1-eta)/(1-eta)+ev(i,iz)
               if(tmpmax&lt;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
  </code></pre>
  
  <p><strong>Matlab parfor:</strong>
  (I know the code below could be made faster by using vectorized syntax, but the whole point of the exercise is to compare loop speeds).</p>
  
  <pre><code>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&gt;tol
  
      parfor ik=1:nk
            for iz = 1:nz
            tmpmax = -intmax;
  
            for i = 1:nk
               c1 = c0(ik,iz) - kgrid(i);
               if (c1&lt;0) 
                   continue
               end 
               c1 = c1^(1-eta)/(1-eta)+ev(i,iz);
               if tmpmax&lt;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
  </code></pre>
  
  <p><strong>Matlab CUDA:</strong></p>
  
  <p>This is what I have no clue how to code. Is using <code>arrayfun</code> the only way of doing this? In fortran is so simple to move from OpenMP to OpenACC. Isn't there an easy way in Matlab of going from parfor to GPUs loops?</p>
  
  <p>The time comparison between codes:</p>
  
  <pre><code>Fortran OpenMP: 83.1 seconds 
  Fortran ACC:    2.4 seconds
  Matlab parfor:  1182 seconds
  </code></pre>
  
  <p>Final remark, I should say the codes above solve a simple Real Business Cycle Model and were written based on <a href="https://www.kellogg.northwestern.edu/rc/workshops/mp2.html" rel="noreferrer">this</a>.</p>