PyCUDA точность кода умножения матрицы

Question

Я пытаюсь изучить CUDA и использовать PyCUDA для написания простого кода умножения матрицы. В течение двух 4x4 матриц генерируются случайным образом я получаю следующее решение:

Cuda: 
[[ -5170.86181641 -21146.49609375 20690.02929688 -35413.9296875 ] 
[-18998.5   -3199.53271484 13364.62890625 7141.36816406] 
[ 31197.43164062 21095.02734375 1750.64453125 11304.63574219] 
[ -896.64978027 18424.33007812 -17135.00390625 7418.28417969]] 

Python: 
[[ -5170.86035156 -21146.49609375 20690.02929688 -35413.9296875 ] 
[-18998.5   -3199.53271484 13364.62695312 7141.36816406] 
[ 31197.43164062 21095.02929688 1750.64404297 11304.63574219] 
[ -896.64941406 18424.33007812 -17135.00390625 7418.28417969]] 

Cuda-Python: 
[[-0.00146484 0.   0.   0.  ] 
[ 0.   0.   0.00195312 0.  ] 
[ 0.   -0.00195312 0.00048828 0.  ] 
[-0.00036621 0.   0.   0.  ]] 

погрешность составляет порядка 1e-3 и возрастает по мере увеличить размер матриц. Я не уверен, что это ошибка или нет. Мой вопрос заключается в том, является ли такая «большая» ошибка нормальной вещью с int32 или я делаю что-то неправильно?

Вот исходный код:

matmul.py

import numpy as np 
import time 
# import pycuda stuff 
import pycuda.driver as cuda 
import pycuda.autoinit 
from pycuda.compiler import SourceModule 

BLOCK_SIZE = 16 

n = 4 
ni = np.int32(n) 

# matrix A 
a = np.random.randn(n, n)*100 
a = a.astype(np.float32) 

# matrix B 
b = np.random.randn(n, n)*100 
b = b.astype(np.float32) 

# matrix B 
c = np.empty([n, n]) 
c = c.astype(np.float32) 

# allocate memory on device 
a_gpu = cuda.mem_alloc(a.nbytes) 
b_gpu = cuda.mem_alloc(b.nbytes) 
c_gpu = cuda.mem_alloc(c.nbytes) 

# copy matrix to memory 
cuda.memcpy_htod(a_gpu, a) 
cuda.memcpy_htod(b_gpu, b) 

# compile kernel 
mod = SourceModule(open("kernels.cu", "r").read()) 

# get function 
matmul = mod.get_function("matmul"); 


# set grid size 
if n%BLOCK_SIZE != 0: 
    grid=(n/BLOCK_SIZE+1,n/BLOCK_SIZE+1,1) 
else: 
    grid=(n/BLOCK_SIZE,n/BLOCK_SIZE,1) 

# call gpu function 
start = time.time() 
matmul(ni, a_gpu, b_gpu, c_gpu, block=(BLOCK_SIZE,BLOCK_SIZE,1), grid=grid); 
end = time.time() 
print "Time: %.5f s"%(end-start) 

# copy back the result 
cuda.memcpy_dtoh(c, c_gpu) 

print np.linalg.norm(c - np.dot(a,b)) 
print c 
print np.dot(a,b) 
print c - np.dot(a,b) 

kernels.cu

__global__ void matmul(int n, const float *A, const float *B, float *C){ 

    int tx = threadIdx.x; 
    int ty = threadIdx.y; 

    int bx = blockIdx.x; 
    int by = blockIdx.y; 

    int row = by*blockDim.y + ty; 
    int col = bx*blockDim.x + tx; 

    if(row < n && col < n){ 
    float val = 0.0; 
    for(int i=0; i<n; ++i){ 
     val += A[row*n + i]*B[n*i + col]; 
    } 
    C[row*n + col] = val; 
    } 
} 

Answer 1

Просто добавить к тому, что Уоррен сказал. Я не думаю, что здесь что-то не так. Вы сравниваете результаты с плавающей запятой, создаваемые двумя разными машинами (CPU и GPU). Они не гарантированно побитообразны для операций на уровне, о котором вы думаете, отчасти потому, что порядок операций на графическом процессоре не обязательно совпадает с порядком операций на графическом процессоре. По мере увеличения размера матриц вы увеличиваете количество значений, суммируемых вместе, и ваша абсолютная ошибка увеличивается, так как вы добавляете кучу небольших битовых ошибок вместе.

В общем, эти соображения всегда должны вступать в игру при сравнении результатов с плавающей точкой. Побитовые идентичные результаты редко можно ожидать от двух разных вычислений. И даже что-то простое, как изменение порядка операций, может сделать его другим вычислением для операций с плавающей запятой. Вы можете прочитать this paper, особенно раздел 2.2.