安装了Anaconda的NumbaPro后出现CUDA断言错误

Question

我正在尝试使用NumbaPro的cuda扩展程序来乘以大型数组矩阵。最终我想要的是将一个大小为NxN的矩阵乘以一个作为1D矩阵提供的对角线矩阵（因此，a.dot(numpy.diagflat(b))相当于a * b）。但是，我遇到了一个没有提供任何信息的断言错误。

只有当我乘以两个1D数组矩阵时，才能避免这种断言错误，但这不是我想做的事情。

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

- brebs

---

你能否检查 _cudadispatch.py 的第191行，看看断言是什么？ - BenC

不幸的是，它只存在编译版本。当我反编译它时，行号是不同的，但从我所看到的来看，harrism是正确的，错误是由于它不是标量输入而引发的。我还想指出，唯一有效的反编译器是在回答这个问题[1]中找到的Decompyle++。此外，文件位置是“〜/anaconda/lib/...”，而不是“/opt/anaconda1anaconda2anaconda3/...” [1]: https://dev59.com/TGsy5IYBdhLWcg3w-i_h - brebs

啊，闭源库的乐趣……那就祝你好运吧！ - BenC

Answer 1

问题在于您正在对一个接受非标量参数的函数使用vectorize。 NumbaPro的想法是，它接受一个标量函数作为输入，并生成一个函数，将标量操作并行应用于向量的所有元素。请参阅NumbaPro文档。

您的函数接受矩阵和向量，这些都不是标量。 [编辑]您可以使用NumbaPro的cuBLAS包装器或编写自己的简单内核函数在GPU上执行您想要的操作。以下是一个演示两者的示例。请注意，您需要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 函数更加智能。它只为每个矩阵元素执行单个乘法，并且永远不会将对角线扩展为完整矩阵。这是我在NVIDIA Tesla K20c GPU上测试N = 1000的结果：

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

- harrism

---

谢谢您指出这一点。您是指cuBLAS下的此页面上的函数吗？我没有看到任何用于简单矩阵乘法的函数。我有什么遗漏吗？http://docs.continuum.io/numbapro/cudalib.html - brebs

一种方法是使用CUBLAS SGEMM，但当您知道一个矩阵是对角矩阵时，矩阵乘法可能过于复杂。另一种方法是编写“CUDA Python”内核--这应该更快，因为您可以针对您的情况进行专门优化。然而，在编写示例时我遇到了一些错误，所以在更新我的答案之前，我正在等待Continuum的帮助。感谢您的接受。 - harrism

没问题。你用的是什么GPU？10000*10000个浮点数是400MB。这两种方法（CUBLAS和自定义内核）中的每一种都会生成3个这样大小的数组。如果你只使用一种方法，你需要1.2GB的空间。如果你同时使用这两种方法（上面的代码没有改变），并且NumbaPro在它们之间没有回收内存，那么它将是这个大小的两倍。 - harrism

唯一增加GPU内存的诀窍就是购买新的GPU。 :) 由于这是一个相当琐碎的计算，您可以将其分成块并迭代地处理这些块。 - harrism

实际上，我的diagproduct内核计算是错误的：它不会将对角向量提升为2D矩阵，因此您只需要2*400MB + 40KB的数据，这应该适合您的1GB GPU。 - harrism

Answer 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

- Remi