GPU Parallelization: Kernel for Value Function Iteration

Dear all,

I'd like to write a matlab function that does Value Function Iteration on GPU.

What I have in mind is very similar to this code written in Julia:

Value Function Iteration under uncertainty

I'm really struggling in writing the kernel for the GPU.

In the Julia-Code, doing so looks like this:

# Write kernel for GPU manually:   
gpu_call(grid, (grid, V, policy, z, P, Float32(alpha), Float32(beta), Float32(delta), Float32(sigma), UInt32(SIZE_GRID), UInt32(SIZE_Z)))  
do state, grid, V, policy, z, P, alpha, beta, delta, sigma, SIZE_GRID, SIZE_Z
# Each kernel executes for one value of the capital grid:
idx = @linearidx grid

What would be the equivalent functions in matlab for

gpu_call( )

and

__ = @linearidx __

The only thing similair to gpu_call that I found was:

KERN = parallel.gpu.CUDAKernel(PTXFILE,CPROTO)   

But this requires code based on C for CUDA or OpenCL, if I understand it correctly. I am unable to deal with such code.

I have the parallel computing toolkit installed.

My very rough sketch of the function (without the part I don't get) looks like this:

    function [V,pol] = VFI_own_gpu_attempt(alpha,beta,delta,eta,z_grid,k_grid,pi_z,tol)
    size_k_grid = size(k_grid,1);
    size_z_grid = length(z_grid);
    
    k_grid_G = gpuArray(k_grid);
    z_grid_G = gpuArray(z_grid);
    pi_z_G = gpuArray(pi_z);
    V0 = ones(size_k_grid,size_z_grid,'gpuArrays');
    V = ones(size_k_grid,size_z_grid,'gpuArrays');
    pol = zeros(size_k_grid,size_z_grid,'gpuArrays');
    
    while abs(V-V0)>tol
    V0 = V;
    % write kernel
    %gpu_call(...)
    
    % each kernel executes for one value of the capital grid
    %idx = @linearidx grid
    
    for i_z = 1:size_z_grid
    F = -Inf;
    pol_i = uint(1)
        for i_k = 1_size_k_grid
        c = z_grid_G(i_z)*k_grid_G(idx)^alpha + (1-delta)*k_grid_G(idx) - k_grid_G(i_k)M
            if c>0
            F0 = ((c)^(1-eta)-1)/(1-eta)
                for j = 1:size_z_grid
                    F1 = F0 + beta*pi_z_G(i_z,j)*V(i_k,j);
                end
            end
            if F1 > F
            F = F1;
            pol_i = uint64(i_k);
            end
        end
    V(idx,i_z) = F;
    pol(idx,i_z) = pol_i;
    
    end

Thanks for any help, tips or advice! :)

Answer 1

[Alessandro]

You can try using the VFI toolkit, which uses gpu parallelization in Matlab.

https://github.com/vfitoolkit