Matlab上用于GPU并行化值函数迭代的写内核

Question

我想编写一个在GPU上执行Value Function Iteration的matlab函数。我的想法与用Julia写的代码非常相似：Value Function Iteration under uncertainty

我真的很努力为GPU编写内核。在Julia-Code中，操作如下所示：

# 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

在matlab中等效的功能是什么gpu_call( )和__ = @linearidx __？

我发现唯一能发现gpu_call的东西是：KERN = parallel.gpu.CUDAKernel(PTXFILE,CPROTO)但这需要基于C的CUDA或OpenCL代码，如果我理解正确的话。我无法处理此类代码。

我已经安装了并行计算工具包。

感谢您的任何帮助，提示或建议！ ：）

编辑：我非常粗略的功能草图（没有我不了解的部分）看起来像这样：

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

编辑2：我用arrayfun尝试了另一种方法。main.m是：

alpha=0.35;
beta= 0.984;
delta=0.01;
eta=2;

tol=10^(-4);

% Not using Tauchen but some transition matrix with 3 grid points for stochastic process
z_grid = [0.4,0.8,1.2];
pi_z = [0.7,0.2,0.1;0.1,0.8,0.1;0.05,0.1,0.85];

% capital grid
k_min = 0.01;
k_max = 20;
n_k = 1000;
k_grid = linspace(k_min,k_max,n_k);

% prepare objects for GPU
k = gpuArray(k_grid);
z = gpuArray(z_grid);
pi_z_G = gpuArray(pi_z);

size_k = size(k,2);
size_z = length(z);

V0 = zeros(size_k,size_z,'gpuArray');
V = ones(size_k,size_z,'gpuArray');
pol = zeros(size_k,size_z,'gpuArray');
diff=max(max(V-V0));

while abs(diff)>tol
V0 = V;
y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

for i_k = 1:size_k
    for i_z = 1:size_z

[V,pol]= arrayfun(@VFI_own_gpu_attempt1,V0,y_k,i_k,i_z, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G);
    end
end
diff=max(max(V-V0));
end
toc

VFI功能为：

function [V,pol] = VFI_own_gpu_attempt1(V0,y_k,i_k,i_z,alpha, beta, delta, eta,size_k,size_z,z,k,pi_z_G)
    low_k=1;
    if(y_k(i_k,i_z) > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k(i_k,i_z), 1);
    end

    if(k(high_k) > y_k(i_k,i_z))
    high_k = high_k -1;
    end

    N_k = high_k; 
%maximization 
F = ((y_k(i_k,i_z)*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G(i_z,:)';
[V(i_k,i_z),pol(i_k,i_z)]=max(F);


end

当尝试使用arrayfun运行该行时，代码停止。错误消息显示：

Error using gpuArray/arrayfun

矩阵尺寸必须一致。

Error in own_gpu_attempt (line 41)
[V,pol]= arrayfun(@VFI_own_gpu_attempt1,V0,y_k,i_k,i_z, alpha, beta, delta, eta,size_k,
size_z,z,k,pi_z_G);

但是，当我将第41行从arrayfun转换为普通的CPU函数时，该函数执行正常。怎么可能？

Answer 1

因为GPU上的arrayfun()要求输入要么是网格（或其特定尺寸），要么是标量。在CPU上时，它将接受矩阵等（例如pi_z_G）作为输入。

您可以在vfitoolkit.com（一个Matlab包）中找到可在GPU上运行的通用实现。它在ValueFnIter_Case1()命令中。

您是使用它还是使用它来了解何时以及如何在GPU上使用arrayfun()，取决于您。

本质上是只是更具体地确定代码的哪些部分实际上是要在网格上评估的函数，而哪些不是。

[编辑：我是vfitoolkit的开发者，但我的意思只是关于arrayfun（）在matlab中的工作方式；添加了此内容，因为这不是我打算发送垃圾邮件。]