Ytukyg

I have been playing around with parallelization both using CUDA 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 arrayfun function. But I might be wrong.

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.

Below, I am sharing three codes and benchmarks:

1. Fortran OpenMP code


2. Fortran CUDA code


3. Matlab parfor code


4. Matlab Cuda (?) this is the one I don't know how to do.

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 CUDA:

Just replace the mainloop syntax on the above code with:

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:
(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).

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:

This is what I have no clue how to code. Is using arrayfun 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?

The time comparison between codes:

Fortran OpenMP: 83.1 seconds 

Fortran CUDA: 2.4 seconds

Matlab parfor: 1182 seconds

Final remark, I should say the codes above solve a simple Real Business Cycle Model and were written based on this.