使用Pycuda块和网格处理大数据

Question

我需要帮助来了解我的块和网格的大小。
我正在构建一个Python应用程序，用于基于scipy进行一些度量计算，比如：欧几里得距离、曼哈顿距离、皮尔逊相关系数、余弦相似度等等。

这个项目叫做 PycudaDistances。

在处理小数组时，它似乎运行得很好。但是当我进行更全面的测试时，不幸的是它没有成功。我下载了movielens数据集（http://www.grouplens.org/node/73）。

使用 Movielens 100k，我声明了一个形状为(943, 1682)的数组。也就是说，用户有943个，评估的电影有1682部。对于那些没有被某个用户评估的电影，我把值设置为0。

但是当数组变得更大时，算法就不再工作。我遇到了以下错误：

pycuda._driver.LogicError: cuFuncSetBlockShape失败：无效值。

在研究这个错误时，我发现Andrew解释说，支持512个线程，如果要处理更大的块，就需要使用块和网格。

我希望能得到一些帮助，让欧几里得距离的算法可以适应从小数组到大数组的处理。

def euclidean_distances(X, Y=None, inverse=True):
    X, Y = check_pairwise_arrays(X,Y)
    rows = X.shape[0]
    cols = Y.shape[0]
    solution = numpy.zeros((rows, cols))
    solution = solution.astype(numpy.float32)

    kernel_code_template = """
    #include <math.h>
    
    __global__ void euclidean(float *x, float *y, float *solution) {

        int idx = threadIdx.x + blockDim.x * blockIdx.x;
        int idy = threadIdx.y + blockDim.y * blockIdx.y;
        
        float result = 0.0;
        
        for(int iter = 0; iter < %(NDIM)s; iter++) {
            
            float x_e = x[%(NDIM)s * idy + iter];
            float y_e = y[%(NDIM)s * idx + iter];
            result += pow((x_e - y_e), 2);
        }
        int pos = idx + %(NCOLS)s * idy;
        solution[pos] = sqrt(result);
    }
    """
    kernel_code = kernel_code_template % {
        'NCOLS': cols,
        'NDIM': X.shape[1]
    }

    mod = SourceModule(kernel_code)

    func = mod.get_function("euclidean")
    func(drv.In(X), drv.In(Y), drv.Out(solution), block=(cols, rows, 1))

    return numpy.divide(1.0, (1.0 + solution)) if inverse else solution

想了解更多细节，请查看：https://github.com/vinigracindo/pycudaDistances/blob/master/distances.py

Answer 1

这个被接受的答案在原则上是对的，但talonmies写的代码有点问题。那一行代码： gdim = ( (dx + (mx>0)) * bdim[0], (dy + (my>0)) * bdim[1]) ) 应该改成： gdim = ( (dx + (mx>0)), (dy + (my>0)) ) 除了多了一个明显的括号外，gdim会产生比你想要的线程多得多。talonmies在他的文字中解释得很对，线程的数量是块大小乘以网格大小。但是他列出的gdim会给你总的线程数，而不是你想要的正确的网格大小。

Answer 2

要设置你的内核执行参数，你需要做两件事（按这个顺序）：

1. 确定块大小

块大小主要受硬件限制和性能的影响。我建议你看看这个回答，里面有更详细的信息。简单来说，你的GPU对每个块能运行的线程总数有个限制，而且它的寄存器、共享内存和本地内存大小也是有限的。你选择的块尺寸必须在这些限制之内，否则内核就无法运行。块大小还会影响内核的性能，你需要找到一个能提供最佳性能的块大小。块大小应该始终是32的整数倍，因为到目前为止，所有支持CUDA的硬件的warp大小都是32。

2. 确定网格大小

对于你展示的这种内核，所需的块数量与输入数据的大小和每个块的尺寸直接相关。

举个例子，如果你的输入数组大小是943x1682，而你选择的块大小是16x16，那么你需要一个59 x 106的网格，这样在内核启动时就会产生944x1696个线程。在这种情况下，输入数据的大小不是块大小的整数倍，你需要修改你的内核，以确保它不会读取越界的数据。一个可能的做法是：

__global__ void euclidean(float *x, float *y, float *solution) {
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    int idy = threadIdx.y + blockDim.y * blockIdx.y;

     if ( ( idx < %(NCOLS)s ) && ( idy < %(NDIM)s ) ) {

        .....
     }
}

启动内核的Python代码可能看起来像这样：

bdim = (16, 16, 1)
dx, mx = divmod(cols, bdim[0])
dy, my = divmod(rows, bdim[1])

gdim = ( (dx + (mx>0)) * bdim[0], (dy + (my>0)) * bdim[1]) )
func(drv.In(X), drv.In(Y), drv.Out(solution), block=bdim, grid=gdim)

这个问题和回答也许能帮助你理解这个过程是如何工作的。

请注意，上述所有代码都是在浏览器中编写的，未经过测试。使用时请自行承担风险。

另外，请注意，这些内容是基于对你代码的非常简短的阅读，可能并不准确，因为你在问题中并没有详细描述代码是如何调用的。