Ошибка Anbaonda's NumbaPro CUDA Assertion Error

Question

Я пытаюсь использовать расширение cuda NumbaPro для умножения массивов больших массивов. В конце концов, я хочу умножить матрицу размера NxN на диагональную матрицу, которая будет подаваться в виде 1D матрицы (таким образом, a.dot (numpy.diagflat (b)), который, как я нашел, является синонимом * 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 

Answer 1

Проблема заключается в использовании vectorize на функцию, которая принимает нескалярных аргументы. Идея с NumbaPro's vectorize заключается в том, что она принимает скалярную функцию как входную и генерирует функцию, которая применяет скалярную операцию параллельно всем элементам вектора. См. NumbaPro documentation.

Ваша функция принимает матрицу и вектор, которые определенно не являются скалярными. [Edit] You может делать то, что вы хотите на графическом процессоре, используя обертку NumbaPro для cuBLAS, или создавая собственную простую функцию ядра. Вот пример, который демонстрирует и то, и другое. Примечание потребует 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 

Время включает в себя все копии и из ГПУ, что является существенным препятствием для небольших матриц. Если мы устанавливаем N 10000 и запустить снова, мы получаем гораздо больше: ускорение

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 

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