PyCUDA矩阵乘法的精度代码

我想学习CUDA并用PyCUDA来编写一个简单的矩阵乘法代码。两个4×4随机生成的矩阵的I得到以下溶液：PyCUDA矩阵乘法的精度代码

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顺序和增加的作为我增加矩阵的大小。我不确定它是否有bug。我的问题是这样一个“大”的错误是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; 
    } 
}

来源

2014-01-15 maverick

*相对*错误的顺序为1e-7，这对于单精度（32位）浮点计算来说是合理的。 –

谢谢！我在精确度的相对和绝对顺序之间感到困惑。 – maverick

只是增加了沃伦说。我不认为这里有什么问题。您正在比较两台不同机器（CPU和GPU）生成的浮点结果。他们不能保证按位在你思考水平的操作相同，部分原因是因为操作的GPU上的顺序是不一定相同的操作对GPU的顺序。当你增加矩阵的大小时，你会增加总计在一起的值的数量，并且你的绝对误差会增加，因为你将一堆小的错误加在一起。

通常，这些类型的考虑因素在比较浮点结果时应始终发挥作用。从两个不同的计算中很难预料到按位相同的结果。即使像改变操作顺序这样简单的事情，也可以使其成为浮点操作的不同计算。您可能想要阅读this paper，特别是第2.2节。

来源

2014-01-15 14:51:14

明白了，谢谢！ – maverick

PyCUDA矩阵乘法的精度代码

回答

相关问题