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

几个通用的cmets：

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

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

1. 对于规范矩阵乘法，我通常希望从结果 (C) 矩阵而不是 A 矩阵中得出计算范围。如果您将自己限制在X*Xt 的情况下，那么我认为您使用A 应该没问题。一般情况下不会。

2. 对我来说很明显您有索引问题。我不会尝试将它们全部整理出来，甚至将它们全部识别出来，但我已经指出了一个问题。由于您的网格大小选择，您的 tid 变量涵盖了 0..1023 的范围，这对于此索引模式不可能是正确的：B[tid, col]（B 的行数等于1024)。

3. 在我看来，多个线程可能会写入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
$