我的pyCuda索引是如何工作的？

Posted 2021-04-21

我正在尝试将一个3d数组加载到pycuda（我要加载图像）。我希望每个线程使用for循环处理单个像素的所有通道（这是算法要求）。

到目前为止，我有这个工作：

from pycuda.compiler import SourceModule
mod = SourceModule(open("./cudacode.cu").read())  

multiply_them = mod.get_function("multiply_them")

rows,cols = 800,400
a = numpy.random.randn(rows,cols,3).astype(numpy.float32)
b = numpy.random.randn(rows,cols,3).astype(numpy.float32)



threads = 20
blocks = 16

dest = numpy.zeros_like(a)
multiply_them(
        drv.Out(dest), drv.In(a), drv.In(b),np.int32(rows), np.int32(cols),  # notice it is rows, cols
        block=(blocks,blocks,1), grid=(rows//blocks,cols//blocks))

print( dest- 2*a*b)
print(np.unique( dest- 2*a*b))

我的cuda代码是：

__global__ void multiply_them(float *dest, float *a, float *b,int cols,int rows) # notice it is cols , rows  and not rows, cols 
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if(row<rows && col <cols){

    dest[3* (row*cols + col)] =   2.0* a[3* (row*cols + col)  ]* b[3* (row*cols + col)  ] ; ## first channel of some pixel
    dest[3* (row*cols + col)+1] = 2.0* a[3* (row*cols + col)+1]* b[3* (row*cols + col)+1] ; ## second channel of that pixel
    dest[3* (row*cols + col)+2] = 2.0* a[3* (row*cols + col)+2]* b[3* (row*cols + col)+2] ; ## 3rd channel of that pixel 


    }
}

请注意，在我的cuda函数中，行和列被切换。这整个代码的工作方式很好。

我可以说我的代码工作正常，因为print(np.unique( dest- 2*a*b))打印0.0

dest[3* (row*cols + col)]索引风格，我从这里发现https://wiki.tiker.net/PyCuda/Examples/GlInterop

我的问题是，为什么需要切换max rows和max cols以使其工作？

什么是在逻辑上更正确的方式去做我想要的？

我是cuda的新手，所以请耐心等待，因为我可能会提出可能非常愚蠢的问题

答案

Numpy数组默认使用row-major ordering。 CUDA块中的线程是numbered，因此threadIdx.x是变化最快的维度，而threadIdx.z是最慢的（功能上等于列主要排序）。因此，如果您希望访问行主要有序数据，就像在默认有序numpy数组中那样，并维护允许memory coalescing的访问顺序，那么您应该“反转”在块中使用x和y维度。如果你的numpy数组是列主要有序（order='F'），那么你就不会这样做。

这就是为什么当你的直觉不应该时你的代码正常工作的原因。