二维 Cuda 网格内核中的 Cupy 索引？答案

【问题标题】：Cupy indexing in 2D Cuda Grid Kernels?二维 Cuda 网格内核中的 Cupy 索引？
【发布时间】：2021-12-05 13:03:55
【问题描述】：

我正在尝试开始使用 Cupy 进行一些 Cuda 编程。我需要编写自己的内核。但是，我正在为 2D 内核苦苦挣扎。似乎 Cupy 没有按我预期的方式工作。下面是一个非常简单的 Numba Cuda 中的 2D 内核示例：

import cupy as cp
from numba import cuda

@cuda.jit
def nb_add_arrs(x1, x2, y):
  i, j = cuda.grid(2)
  if i < y.shape[0] and j < y.shape[1]:
    y[i, j] = x1[i, j] + x2[i, j]

x1 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
x2 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
y = cp.zeros((5, 5), dtype=cp.int32)
# Grid and block sizes
tpb = (16, 16)
bpg = (x1.shape[0] // tpb[0] + 1, x1.shape[1] // tpb[0] + 1)
# Call kernel
nb_add_arrs[bpg, tpb](x1, x2, y)

结果如预期：

y
[[2 2 2 2 2]
 [2 2 2 2 2]
 [2 2 2 2 2]
 [2 2 2 2 2]
 [2 2 2 2 2]]

但是，当我尝试在 Cupy 中做这个简单的内核时，我没有得到相同的结果。

cp_add_arrs = cp.RawKernel(r'''
extern "C" __global__
void add_arrs(const float* x1, const float* x2, float* y, int N){
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  int j = blockDim.y * blockIdx.y + threadIdx.y;

  if(i < N && j < N){
    y[i, j] = x1[i, j] + x2[i, j];
  }
}
''', 'add_arrs')

x1 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
x2 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
y = cp.zeros((5, 5), dtype=cp.int32)
N = x1.shape[0]
# Grid and block sizes
tpb = (16, 16)
bpg = (x1.shape[0] // tpb[0] + 1, x1.shape[1] // tpb[0] + 1)
# Call kernel
cp_add_arrs(bpg, tpb, (x1, x2, y, cp.int32(N)))

y
[[0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]

谁能帮我找出原因？

【问题讨论】：

您的 C++ 内核已损坏。如果 cupy 向您显示编译器警告，请不要忽略它们。这：y[i, j] = x1[i, j] + x2[i, j]; 没有做你认为的事情
@talonmies 感谢您的回复。你知道 de C++ 内核有什么问题吗？我遵循了在 numba 内核中使用的相同逻辑。 Cupy 没有显示任何警告。
您需要学习 C++ 才能为 cupy 编写 C++ 内核代码。这种类型的语法[i, j] 通常在 python 中用于二维索引。它在 C++ 中不起作用。你需要做类似[i][j]的事情。
是的，我刚刚告诉你出了什么问题。您不能在 C++ 代码中使用 python 样式的数组索引。正如您已经发现的那样，它什么也不做
请添加您的解决方案作为问题的答案

标签： python cuda cupy

【解决方案1】：

C 中的内存以行优先顺序存储。所以，我们需要按照这个顺序进行索引。此外，由于我传递的是 int 数组，因此我更改了内核的参数类型。代码如下：

cp_add_arrs = cp.RawKernel(r'''
extern "C" __global__
void add_arrs(int* x1, int* x2, int* y, int N){
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  int j = blockDim.y * blockIdx.y + threadIdx.y;
  
  if(i < N && j < N){
    y[j + i*N] = x1[j + i*N] + x2[j + i*N];
  }
}
''', 'add_arrs')

x1 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
x2 = cp.ones(25, dtype=cp.int32).reshape(5, 5)
y = cp.zeros((5, 5), dtype=cp.int32)
N = x1.shape[0]
# Grid and block sizes
tpb = (16, 16)
bpg = (x1.shape[0] // tpb[0] + 1, x1.shape[1] // tpb[0] + 1)
# Call kernel
cp_add_arrs(bpg, tpb, (x1, x2, y, cp.int32(N)))

y
array([[2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2]], dtype=int32)

【讨论】：