PyCUDA核函数

Question

我是PyCUDA的新手，在PyCUDA网站上浏览了一些例子。我正试图找出某些代码行背后的逻辑，如果有人解释其背后的想法，我将非常感激。

下面的代码片段来自PyCUDA网站。在函数定义中，我不明白

int idx = threadIdx.x + threadIdx.y*4;

如何使用上面的行计算数组的索引。为什么会将线程Idx.x和线程Idx.y相加在一起，以及为什么线程Idx.y乘以4。

对于对GPU的函数调用，为什么将块定义为5,5,1。因为它是一个由5x5元素组成的数组，所以在我的理解中，块大小应该是5,5，而不是5,5,1。

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
a = numpy.random.randn(5,5)
a = a.astype(numpy.float32)
a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)
mod = SourceModule("""
__global__ void doubleMatrix(float *a)
{
int idx = threadIdx.x + threadIdx.y*4;
a[idx] *= 2;
}
""")
func = mod.get_function("doubleMatrix")
func(a_gpu, block=(5,5,1))
a_doubled = numpy.empty_like(a)
cuda.memcpy_dtoh(a_doubled, a_gpu)
print ("ORIGINAL MATRIX")
print a
print ("DOUBLED MATRIX AFTER PyCUDA EXECUTION")
print a_doubled

Answer 1

您发布的示例似乎来自(或被抄袭到)一本名为"Python并行编程Cookbook“的书中，我直到五分钟前才听说过它。老实说，如果我是那本书的作者，我会感到羞愧，因为我包含了这样一个陈词滥调的、破碎的例子。

下面是对您发布的内容及其输出的一个小修改：

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
a = numpy.random.randn(5,5)
a = a.astype(numpy.float32)
a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)
mod = SourceModule("""
__global__ void doubleMatrix(float *a)
{
     int idx = threadIdx.x + threadIdx.y*4;
     a[idx] *= 2.f;
}
""")
func = mod.get_function("doubleMatrix")
func(a_gpu, block=(5,5,1))
a_doubled = numpy.empty_like(a)
cuda.memcpy_dtoh(a_doubled, a_gpu)
print a_doubled - 2.0*a

警告: Python 2语法

In [2]: %run matdouble.py
[[ 0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.        ]
 [ 0.         -0.62060976  0.49836278 -1.60820103  1.71903515]]

也就是说，代码不像预期的那样工作，这可能是造成混乱的根源。

this very recent answer中描述了解决存储在线性内存中的多维数组(如numpy数组)的正确方法。任何明智的程序员都会在示例中编写内核，如下所示：

__global__ void doubleMatrix(float *a, int lda)
{
     int idx = threadIdx.x + threadIdx.y * lda;
     a[idx] *= 2.f;
}

这样数组的前导维度作为参数传递给内核(在本例中应该是5，而不是4)。这样做会产生以下结果：

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
a = numpy.random.randn(5,5)
a = a.astype(numpy.float32)
a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)
mod = SourceModule("""
__global__ void doubleMatrix(float *a, int lda)
{
     int idx = threadIdx.x + threadIdx.y * lda;
     a[idx] *= 2.f;
}
""")
func = mod.get_function("doubleMatrix")
lda = numpy.int32(a.shape[-1])
func(a_gpu, lda, block=(5,5,1))
a_doubled = numpy.empty_like(a)
cuda.memcpy_dtoh(a_doubled, a_gpu)
print a_doubled - 2.0*a

产生预期结果：

In [3]: %run matdouble.py
[[ 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.]]

对于对GPU的函数调用，为什么将块定义为5,5,1。因为它是一个由5x5元素组成的数组，所以在我的理解中，块大小应该是5,5，而不是5,5,1。

在数据自动化系统( CUDA )中，所有区块都隐含着三个维度。块大小(5,5)与块大小(5,5,1)相同。最后一个维度可以忽略，因为它是一个维度(即块中的所有线程都有threadIdx.z = 1)。您不应该陷入这样的陷阱:将CUDA块或网格的维度与输入数组的维度混为一谈。有时候让它们保持不变是很方便的，但同样地，这样做也是不必要的，甚至是不可取的。本例中正确编写的BLAS风格内核(假设主要存储顺序行)可能如下所示：

__global__ void doubleMatrix(float *a, int m, int n, int lda)
{
     int col = threadIdx.x + blockIdx.x * blockDim.x;
     int row = threadIdx.y + blockDim.y * blockDim.y;

     for(; row < m; row += blockDim.y * gridDim.y) {
         for(; col < n; col += blockDim.x * gridDim.x) {
             int idx = col + row * lda;
             a[idx] *= 2.f;
         }
    }
}

注:用浏览器编写，未编译或测试

在这里，任何合法的块和网格维度都将正确地处理任意大小的输入数组，元素总数将被放入一个带符号的32位整数中。如果运行过多的线程，有些线程将什么也不做。如果运行的线程太少，一些线程将处理多个数组元素。如果您运行的网格与输入数组具有相同的维度，那么每个线程将处理一个输入，就像您正在学习的示例中的意图一样。如果您想了解如何选择最合适的块和网格大小，我建议启动here。