Anaconda 的 NumbaPro CUDA 断言错误答案

【问题标题】：Anaconda's NumbaPro CUDA Assertion ErrorAnaconda 的 NumbaPro CUDA 断言错误
【发布时间】：2013-06-13 18:55:31
【问题描述】：

我正在尝试使用 NumbaPro 的 cuda 扩展来乘以大型数组矩阵。我最终想要的是将大小为 NxN 的矩阵乘以一个对角矩阵，该矩阵将作为一维矩阵输入（因此，a.dot(numpy.diagflat(b)) 我发现它是 a 的同义词* b)。但是，我收到了一个未提供任何信息的断言错误。

如果我将两个一维数组矩阵相乘，我只能避免这个断言错误，但这不是我想要做的。

from numbapro import vectorize, cuda
from numba import f4,f8
import numpy as np

def generate_input(n):
    import numpy as np
    A = np.array(np.random.sample((n,n)))
    B = np.array(np.random.sample(n) + 10)
    return A, B

def product(a, b):
    return a * b

def main():
    cu_product = vectorize([f4(f4, f4), f8(f8, f8)], target='gpu')(product)

    N = 1000

    A, B = generate_input(N)
    D = np.empty(A.shape)

    stream = cuda.stream()

    with stream.auto_synchronize():
        dA = cuda.to_device(A, stream)
        dB = cuda.to_device(B, stream)
        dD = cuda.to_device(D, stream, copy=False)
        cu_product(dA, dB, out=dD, stream=stream)
        dD.to_host(stream)

if __name__ == '__main__':
    main()

这是我的终端输出的内容：

Traceback (most recent call last):
  File "cuda_vectorize.py", line 32, in <module>
    main()
  File "cuda_vectorize.py", line 28, in main
    cu_product(dA, dB, out=dD, stream=stream)
  File "/opt/anaconda1anaconda2anaconda3/lib/python2.7/site-packages/numbapro/_cudadispatch.py", line 109, in __call__
  File "/opt/anaconda1anaconda2anaconda3/lib/python2.7/site-packages/numbapro/_cudadispatch.py", line 191, in _arguments_requirement
AssertionError

【问题讨论】：

你不能在第 191 行检查 _cudadispatch.py 以了解断言的确切内容吗？
不幸的是它只作为编译版本存在。当我未编译它时，行号是不同的，但据我所知，harrism 是正确的，因为它不是标量输入而导致错误被抛出。我还想指出唯一有效的反编译器是在回答这个问题时发现的 Decompyle++ [1]。此外，文件位置是“~/anaconda/lib/...”而不是“/opt/anaconda1anaconda2anaconda3/...”[1]：stackoverflow.com/questions/8189352/decompile-python-2-7-pyc
啊，封闭源代码库的乐趣...祝你好运！

标签： python numpy anaconda numba numba-pro

【解决方案1】：

问题是你在一个接受非标量参数的函数上使用vectorize。 NumbaPro 的vectorize 的想法是，它将标量函数作为输入，并生成一个函数，该函数将标量运算并行应用于向量的所有元素。请参阅NumbaPro documentation。

你的函数需要一个矩阵和一个向量，它们绝对不是标量。 [编辑] 您可以在 GPU 上使用 NumbaPro 的 cuBLAS 包装器或编写您自己的简单内核函数来做您想做的事情。这是一个演示两者的示例。 Note 需要 NumbaPro 0.12.2 或更高版本（本次编辑刚刚发布）。

from numbapro import jit, cuda
from numba import float32
import numbapro.cudalib.cublas as cublas
import numpy as np
from timeit import default_timer as timer

def generate_input(n):
    A = np.array(np.random.sample((n,n)), dtype=np.float32)
    B = np.array(np.random.sample(n), dtype=A.dtype)
    return A, B

@cuda.jit(argtypes=[float32[:,:], float32[:,:], float32[:]])
def diagproduct(c, a, b):
  startX, startY = cuda.grid(2)
  gridX = cuda.gridDim.x * cuda.blockDim.x;
  gridY = cuda.gridDim.y * cuda.blockDim.y;
  height, width = c.shape

  for y in range(startY, height, gridY):
    for x in range(startX, width, gridX):       
      c[y, x] = a[y, x] * b[x]

def main():

    N = 1000

    A, B = generate_input(N)
    D = np.empty(A.shape, dtype=A.dtype)
    E = np.zeros(A.shape, dtype=A.dtype)
    F = np.empty(A.shape, dtype=A.dtype)

    start = timer()
    E = np.dot(A, np.diag(B))
    numpy_time = timer() - start

    blas = cublas.api.Blas()

    start = timer()
    blas.gemm('N', 'N', N, N, N, 1.0, np.diag(B), A, 0.0, D)
    cublas_time = timer() - start

    diff = np.abs(D-E)
    print("Maximum CUBLAS error %f" % np.max(diff))

    blockdim = (32, 8)
    griddim  = (16, 16)

    start = timer()
    dA = cuda.to_device(A)
    dB = cuda.to_device(B)
    dF = cuda.to_device(F, copy=False)
    diagproduct[griddim, blockdim](dF, dA, dB)
    dF.to_host()
    cuda_time = timer() - start   

    diff = np.abs(F-E)
    print("Maximum CUDA error %f" % np.max(diff))

    print("Numpy took    %f seconds" % numpy_time)
    print("CUBLAS took   %f seconds, %0.2fx speedup" % (cublas_time, numpy_time / cublas_time)) 
    print("CUDA JIT took %f seconds, %0.2fx speedup" % (cuda_time, numpy_time / cuda_time))

if __name__ == '__main__':
    main()

内核明显更快，因为 SGEMM 执行完整的矩阵-矩阵乘法 (O(n^3))，并将对角线扩展为完整矩阵。 diagproduct 函数更智能。它只是对每个矩阵元素进行一次乘法运算，并且从不将对角线扩展为完整矩阵。以下是我在 N=1000 的 NVIDIA Tesla K20c GPU 上的结果：

Maximum CUBLAS error 0.000000
Maximum CUDA error 0.000000
Numpy took    0.024535 seconds
CUBLAS took   0.010345 seconds, 2.37x speedup
CUDA JIT took 0.004857 seconds, 5.05x speedup

时间包括进出 GPU 的所有副本，这对于小型矩阵来说是一个重大瓶颈。如果我们将 N 设置为 10,000 并再次运行，我们将获得更大的加速：

Maximum CUBLAS error 0.000000
Maximum CUDA error 0.000000
Numpy took    7.245677 seconds
CUBLAS took   1.371524 seconds, 5.28x speedup
CUDA JIT took 0.264598 seconds, 27.38x speedup

但是，对于非常小的矩阵，CUBLAS SGEMM 具有优化的路径，因此更接近 CUDA 性能。这里，N=100

Maximum CUBLAS error 0.000000
Maximum CUDA error 0.000000
Numpy took    0.006876 seconds
CUBLAS took   0.001425 seconds, 4.83x speedup
CUDA JIT took 0.001313 seconds, 5.24x speedup

【讨论】：

感谢您指出这一点。你指的是cuBLAS下这个页面上的功能吗？我没有看到任何用于简单矩阵乘法的函数。我错过了什么吗？ docs.continuum.io/numbapro/cudalib.html
一种方法是 CUBLAS SGEMM，但是当您知道一个矩阵是对角线时，矩阵-矩阵乘法就显得多余了。另一个是编写“CUDA Python”内核——这应该更快，因为您可以针对您的情况专门使用它。但是，我在编写示例时遇到了错误，因此在更新答案之前，我正在等待 Continuum 的帮助。感谢您的接受。
好的，我用一个工作示例更新了我的答案。感谢 Continuum Analytics 快速修复了我在此过程中发现的错误！
非常感谢！这不仅在提高速度方面非常有用，而且对于学习如何将 cuda 应用到我的函数中也非常有用。不幸的是，我在 N=8000 之后得到了一个 CUDA_ERROR_OUT_OF_MEMORY。我的目标是在 N=15,000 的矩阵上运行它。你知道这需要多少内存/运行 N=10,000 需要多少内存？
不用担心。你有什么GPU？ 10000*10000 浮点数为 400MB。这两种方法（CUBLAS 和自定义内核）中的每一种都会生成 3 个该大小的数组。如果只使用一种方法，则需要 1.2GB。如果您同时使用这两种方法（上面的代码，未更改），并且如果 NumbaPro 不回收它们之间的内存，那将是两倍。

【解决方案2】：

只是为了重新考虑所有这些考虑因素。我还想在 CUDA 上实现一些矩阵计算，但后来听说了 numpy.einsum 函数。事实证明，einsum 的速度非常快。在这种情况下，这里是它的代码。但它可以应用于多种类型的计算。

G = np.einsum('ij,j -> ij',A, B)

在速度方面，这里是 N = 10000 的结果

Numpy took    8.387756 seconds
CUDA JIT took 0.218394 seconds, 38.41x speedup
EINSUM took 0.131751 seconds, 63.66x speedup

【讨论】：