cuda.jit 矩阵乘法崩溃

Question

我正在尝试编写 cuda.jit 矩阵乘法，并限制我的线程块数量，它只能是一个。而且我也知道我的乘法是 X*Xtranspose 的形式。

    def matmul_gpu(X, Y):
    # Allocate the output matrix in GPU memory using cuda.to_device
    #
    # invoke the dot kernel with 1 threadBlock with 1024 threads
    #
    # copy the output matrix from GPU to cpu using copy_to_host()
    gpu_mat1 = cuda.to_device(X)
    gpu_mat2 = cuda.to_device(Y)
    res = np.zeros(shape=(X.shape[0], Y.shape[1]), dtype=np.float32)
    gpu_mult_res = cuda.to_device(res)
    threads_per_block = 1024
    blocks_per_grid = 1
    matmul_kernel[blocks_per_grid, threads_per_block](
        gpu_mat1, gpu_mat2, gpu_mult_res)
    mult_res = gpu_mult_res.copy_to_host()
    return mult_res


@cuda.jit
def matmul_kernel(A, B, C):
    num_of_threads = cuda.gridsize(1)
    tid = cuda.grid(1)
    rows_num = A.shape[0]
    cols_num = A.shape[1]
    step = int(np.math.ceil(num_of_threads / cols_num))
    row = int(np.math.floor(tid / cols_num))
    col = int(tid % cols_num)
    for row_start_idx in range(0, rows_num, step):
        if row_start_idx + row < rows_num and col < cols_num:
            C[row_start_idx + row, col] += A[row_start_idx + row, tid] * B[tid, col]

对于维度为 128,256 或 256,128 的矩阵，它会崩溃，并且会按此顺序通过回溯抛出这些错误。

...
 Call to cuMemcpyDtoH results in UNKNOWN_CUDA_ERROR
...
Call to cuMemFree results in UNKNOWN_CUDA_ERROR

它适用于非常大的尺寸，如 1024、2048 和 2048、1024，并且它适用于具有相同尺寸的输入，但有时尺寸不同会引发上述错误。除了我刚刚注意到的 256*256 之外，它几乎从不抛出相等尺寸的错误，所以它应该与那些相关。

调试帮助代码：

# this is the comparison function - keep it as it is, don't change X or Y.
def matmul_comparison():
    X = np.random.randn(1000, 1024)
    Y = np.random.randn(1024, 1000)

    def timer(f):
        return min(timeit.Timer(lambda: f(X, Y)).repeat(3, 5))

    # print('Python:', timer(matmul_trivial)) we will not consider this since it takes infinite time :)
    #print('Numpy:', timer(np.matmul))
    #print('Numba:', timer(matmul_numba))
    print('CUDA:', timer(matmul_gpu))


if __name__ == '__main__':
    os.environ['NUMBAPRO_NVVM'] = '/usr/local/cuda-9.0/nvvm/lib64/libnvvm.so'
    os.environ['NUMBAPRO_LIBDEVICE'] = '/usr/local/cuda-9.0/nvvm/libdevice/'
    matmul_comparison()

Answer 1

一些一般性评论：

除非我验证了数字的正确性，否则我不会宣布一切正常。
我认为进行错误检查是一种很好的做法。 cuda-memcheck 工具可以检查各种错误，即使使用 numba python CUDA 代码也是如此。您的代码会抛出错误，即使您建议的尺寸可以正常工作。
朴素的矩阵乘法具有相当典型的格式，并且在很多地方都有介绍，例如 the CUDA programming guide。如果我正在研究这个，如果可能的话，我会以此为起点。
从性能的角度来看，在 1024 个线程的单个线程块上任意限制 CUDA 代码运行是一个非常糟糕的主意。我无法想象你为什么要那样做。
尽管如此，如果我想使用任意网格排列来处理 CUDA 算法，规范技术将是 grid-stride loop。

关于您的代码，一些问题立即出现：

对于规范矩阵乘法，我通常希望从结果 (C) 矩阵而不是 A 矩阵中得出计算范围。如果您将自己限制在 X*Xt 情况下，那么我认为您使用 A 应该没问题。一般情况下不是。
对我来说很明显你有索引问题。我不会尝试将它们全部分类或什至全部识别出来，但我已经指出了一个问题。由于您选择的网格大小，您的 tid 变量涵盖了 0..1023 的范围，并且对于此索引模式来说这不可能是正确的：B[tid, col]（除了 [= 的行数18=] 等于 1024)。
在我看来，您有可能让多个线程写入 C 矩阵中的同一输出位置。 CUDA 不会为您解决这个问题。您不应该期望让多个线程写入同一输出位置才能正常工作，除非您已采取措施通过原子或经典并行缩减来实现这一点。而且我不想将这些方法中的任何一种引入到这个问题中，所以我认为基本方法很麻烦。

可能还有其他问题。但是由于上面的考虑 3，与其尝试修复您的代码，我宁愿从规范的朴素矩阵乘法开始并使用网格步幅循环。

这是一个结合了这些想法的例子：

$ cat t59.py
import numpy as np
from numba import cuda,jit


@cuda.jit
def matmul_kernel(A, B, C):
    num_of_threads = cuda.gridsize(1)
    tid = cuda.grid(1)
    rows_num = C.shape[0]
    cols_num = C.shape[1]
    idx_range = A.shape[1]
    for mid in range(tid, rows_num*cols_num, num_of_threads):
        row = mid // cols_num
        col = mid - (row*cols_num)
        my_sum = 0.0
        for idx in range(0, idx_range):
            my_sum += A[row, idx] * B[idx, col]
        C[row, col] = my_sum

def matmul_gpu(X, Y):
    # Allocate the output matrix in GPU memory using cuda.to_device
    #
    # invoke the dot kernel with 1 threadBlock with 1024 threads
    #
    # copy the output matrix from GPU to cpu using copy_to_host()
    gpu_mat1 = cuda.to_device(X)
    gpu_mat2 = cuda.to_device(Y)
    res = np.zeros(shape=(X.shape[0], Y.shape[1]), dtype=np.float32)
    gpu_mult_res = cuda.to_device(res)
    threads_per_block = 1024
    blocks_per_grid = 1
    matmul_kernel[blocks_per_grid, threads_per_block](
        gpu_mat1, gpu_mat2, gpu_mult_res)
    mult_res = gpu_mult_res.copy_to_host()
    return mult_res

wA = 256
hA = 128
wB = hA
hB = wA


mA = np.ones(shape=(hA,wA), dtype=np.float32)
mB = np.ones(shape=(hB,wB), dtype=np.float32)
mC = matmul_gpu(mA,mB)
print(mC)
$ cuda-memcheck python t59.py
========= CUDA-MEMCHECK
[[ 256.  256.  256. ...,  256.  256.  256.]
 [ 256.  256.  256. ...,  256.  256.  256.]
 [ 256.  256.  256. ...,  256.  256.  256.]
 ...,
 [ 256.  256.  256. ...,  256.  256.  256.]
 [ 256.  256.  256. ...,  256.  256.  256.]
 [ 256.  256.  256. ...,  256.  256.  256.]]
========= ERROR SUMMARY: 0 errors
$

cuda.jit 矩阵乘法崩溃

cuda.jit matrices multiplication crashes

jit

cuda

gpu

matrix-multiplication

numba