Solve RBC Model with GPU computing?

Dear all

I'm trying to solve a simple RBC model, similar to the one in Aldrich (2010), with GPU computing.

In the paper, they solve it in OpenCL and CUDA, but since then Matlab added some GPU functionalities.

I'm now trying to figure out, how difficult it is to speed up the value function iteration inside Matlab, without using any OpenCl or CUDA.

These points, I don't quite understand at the moment:

• if I know how to implement the VFI for the RBC model with basic matlab commands, how much additional knowledge do I need to implement it with GPU compute? How much do I have to change the logic of the value function iteration, to adapt to GPU computing?
• -do I have to write a kernel, or can I do the VFI with the in-built functions (e.g. gpuArrays)?

I tried to do this by putting different vectors/matrizes in the

a_g =gpuArray(a)

form. This worked but slowed down the VFI a lot.

EDIT:

Here's the code I' d like to improve upon: [The most important part is the iteration in the last block. Especially the r(row,i)=... line is computationally demanding. I know that I can improve the point with vectorization but would like to try gpu computing on this relatively easy example first.]

%% u(c)=ln(c) ; y=f(k)=k^alpha 

alpha= 0.33; beta= 0.96; %parameter
mink= 0.001; maxk= 0.3 ;% bounds for capital grid
nk=500; % number of grid points

% steady state
kss=(1/(alpha*beta))^(1/(alpha-1));
css=kss^alpha-kss;

% iterate on value function
kgrid = linspace(mink,maxk,nk); % create capital grid
v0= zeros(nk,1); %first guess
v1= zeros(nk,1); %updated guess
vij= zeros(nk,1); %needed for maximization
tol=1e-4; %error tolerance
ii=0;
maxit=100;
%vss=log(f(kss)-kss)+beta*vss, hence vss=log(f(kss)-kss)/(1-beta) is initial guess
cssvector= repmat( css, [nk,1] ); %needed for first guess
v0=log(cssvector)/(1-beta); %first guess
v1= zeros(nk,1); %updated guess
tol=1e-3; %error tolerance
ii=0;
error=inf;

R=zeros(nk,nk); 

for i=1:nk,     % calculate R; attention: for every possible k, theres a column of possible c's
   
     R(:,i)=kgrid(i)^alpha-kgrid;%choices of c    
end
 R(R< 0) =R(R<0)-99999999999; %make negative values super negative
 R(R>0)=log(R(R>0)); %utilities of positive values
 row=[];
 r=R; %needed for iteration 
 
 tic
while (ii<=maxit)
   for i=1:nk,  %value function for all k, only iterate on R where consumption was not negative
    row=find(r(:,i)>-10000);
    r(row,i)=r(row,i)+beta*v0(row);
    v1(i) = max(r(:,i));  % find best choice for K_j
   end
   r=R;
   v0=v1;
  ii=ii+1  
end
toc;

Thanks a lot for any help! :)