矩阵乘法代码的PyCUDA精度答案

【问题标题】：PyCUDA precision of matrix multiplication code矩阵乘法代码的PyCUDA精度
【发布时间】：2014-01-15 22:51:14
【问题描述】：

我正在尝试学习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;
  }
}

【问题讨论】：

relative 误差为 1e-7 级，对于单精度（32 位）浮点计算是合理的。
谢谢！我对准确度的相对顺序和绝对顺序感到困惑。
@maverick ：我正在学习 cuda 编程，你的例子看起来很整齐。我尝试运行代码，但出现错误。 “TypeError：没有注册的转换器能够从这个浮点类型的 Python 对象生成 unsigned int 类型的 C++ 右值”

标签： python cuda pycuda

【解决方案1】：

只是补充沃伦所说的话。我不认为这里有什么问题。您正在比较两台不同机器（CPU 和 GPU）生成的浮点结果。对于您正在考虑的级别的操作，它们不能保证按位相同，部分原因是 GPU 上的操作顺序不一定与 GPU 上的操作顺序相同。随着矩阵大小的增加，加在一起的值的数量也在增加，并且绝对误差也会增加，因为您将一堆小的位错误加在一起。

通常，在比较浮点结果时，应始终考虑这些类型的考虑因素。从两个不同的计算中很少会期望按位相同的结果。对于浮点运算，即使是像改变运算顺序这样简单的事情也可以使它成为不同的计算。您可能想阅读this paper，尤其是第 2.2 节。

【讨论】：