Pycuda中的块和网格如何处理大数据

Question

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

该项目是PycudaDistances。

它似乎在小数组中运行得非常好。当我执行更详尽的测试时，不幸的是它没有工作。我下载了movielens集(http://www.grouplens.org/node/73)。

使用Movielens 100k，我声明了一个形状为(943,1682)的数组。也就是说，用户有943个，评估了1682部电影。未被分类器用户评估的电影我将其配置为0。

使用更大的数组算法不再起作用。我面临以下错误: pycuda._driver.LogicError: cuFuncSetBlockShape failed: invalid value.

研究这个错误，我发现需要使用块和网格来处理更大的块。

我希望你能帮忙将欧几里得距离算法适应于从小到大的数组。

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

- Vinnicyus Gracindo

Answer 1

为了确定你的内核执行参数大小，你需要按照以下两个步骤进行（顺序不能颠倒）：

1. 确定块大小

块大小很大程度上由硬件限制和性能决定。我建议阅读这篇答案以获取更详细的信息，但简短的总结是：你的GPU对每个块可以运行的总线程数有一个限制，并且它有一个有限的寄存器文件、共享和本地内存大小。你选择的块维度必须在这些限制范围内，否则内核将无法运行。块大小也会影响内核的性能，你会发现某个块大小可以给出最佳性能。块大小应该始终是warp大小（在所有支持CUDA的硬件中都是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)

这个问题和答案也可以帮助理解这个过程是如何工作的。

请注意，上面所有的代码都是在浏览器中编写的，并且从未经过测试。使用它时需自担风险。

还请注意，这是基于对您的代码进行了非常简短的阅读，并且可能不正确，因为您在问题中并没有真正描述代码是如何被调用的。

- talonmies

Answer 2

接受的答案在原则上是正确的，然而talonmies列出的代码不完全正确。这行：gdim = ((dx + (mx>0)) * bdim[0], (dy + (my>0)) * bdim[1])应该改为：gdim = ((dx + (mx>0)), (dy + (my>0)))除了一个明显的额外括号外，gdim会产生比你想要的线程数多得多的线程。talonmies在他的文本中已经解释了线程是块大小*网格大小。然而他列出的gdim会给你总线程数而不是所需的正确网格大小。

- DLH