在 Matlab 上为 GPU 并行化值函数迭代编写内核

Question

我想编写一个在 GPU 上执行值函数迭代的 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 中的工作原理；已添加此内容，因为我无意发送垃圾邮件。]