使用特定输入的cuda/cublas简单内核中的数值错误

Numerical error in cuda/cublas simple kernel using particular input

本文关键字：内核 简单 错误 cublas 输入 cuda      更新时间：2023-10-16

我正在使用cuda和cublas，我试图实现简单的运算，如矩阵元素乘法/除法。我的实验只使用浮子。我知道最明显的方法是写一个像这样的内核：

__global__ void mul_elementwise(const unsigned int n, float* source, float* dest, const float value)
{
const unsigned int offset = blockIdx.x * blockDim.x + threadIdx.x;
const unsigned int stride = blockDim.x * gridDim.x;
for (unsigned int i = offset; i < n; i += stride)
{
dest[i] = source[i] * value;
}
}

这个内核可以同时用于乘法和除法(只需使用1/x作为值)。但这也可以使用cublas库来实现：假设我们有一个以列主样式存储的矩阵a m x n和一个标量x，然后将alpha=x或alpha=1/x和d_ones设置为m*n 1s的向量，我们可以调用并获得相同的结果

cublasSaxpy(cublas_handle, m * n, &alpha, d_ones, 1, A_dev, 1);

这两种方法都可以很好地工作，但我在某些特定的矩阵中遇到了一些问题，两种方法对此都不起作用。我隔离了这个大矩阵，并在这里构建了一个MCVE(你可以用nvcc mcve.cu -lcublas编译它。正如你所看到的，这两种情况下的结果都是完全错误的：主机结果完全不同，我正在努力弄清楚发生了什么。我在代码中没有看到任何错误，但也许我应该尝试使用double而不是float，看看会发生什么。对这种情况有什么看法吗？提前感谢！

EDIT#1我尝试过使用doubles，但如果我使用cublasDaxpy，则没有任何变化，同时它与自定义内核完美配合。我认为这些值太小，所以单浮点精度是不够的。

有趣的MCVE。难道不可能将向量缩小到只有几个元素吗？难道不能仅基于一个矢量元素来显示计算差异吗？

不管怎样，我看到了几个问题。

  1. 您的内核实现了以下函数：y=alpha*x。但是SAXPY实现了y=alpha*x+y。现在，如果y一开始是(全部)零，那么这两个将是相同的。但这不是你所拥有的：

    CUBLAS    Your Kernel
---------------------------
alpha:   alpha     alpha
x:       1     ahost    (ahost is  your huge data array)
y:   ahost         -

    所以您的内核正在计算y=alpha * ahost，但是您的CUBLAS调用正在计算y = alpha*1 + ahost。总的来说，我不会期望从这些方面得到同样的结果。

  2. 你对错误的分析在某些方面似乎有缺陷。首先，您要计算float变量中的绝对误差(一个始终为正的数字，因为它是绝对值)，然后将其与负数字进行比较：

    float diff = abs(host[i]-dev[i]);
...
if (diff > (-1e12))

    if测试不总是真的吗？也许你指的是1e-12，尽管这仍然有缺陷。在浮点比较中寻找固定误差阈值应按被比较数字的大小进行缩放。float数量只包含大约6-7个精确的十进制数字。(对这些误差求和也很麻烦。)

这是一个完整的代码，修复了上述问题，并对所有比较(主机<->kernel和主机<->cublas)产生零和错误：

static float array[] = {0x00000000,
0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0xB58DA1CF,0xB50D2FEC,0x34A48536,0xB4A1D5BC,0x358E1345,0x35943AAC,0xB5983F40,0xB43628BB,0xB4A95348,0xB4DB751C,0xB50C8D1A,0xB3EFCBB5,0x3552B8CD,0x3538A167,0x358FDE0D,0xB4D54CE9,0xB5D29BB7,0xB4A234EE,0x346EF2F4,0x35B5D9F2,0xB40F1487,0x3554BC20,0x33FD9466,0xB536D37D,0xB3C2E594,0xB59DA581,0x3584FC87,0x34438F09,0x35D293CB,0xB4FBB002,0xB59F41E9};
#include <iostream>
#include <stdio.h>
#include <cublas_v2.h>
#include <assert.h>
#define TOL 0.0001
typedef unsigned int u32;
#define GET_STRIDE() u32(blockDim.x * gridDim.x)
#define GET_OFFSET() u32(blockIdx.x * blockDim.x + threadIdx.x)
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
if (result != cudaSuccess) {
fprintf(stderr, "CUDA Runtime Error: %sn", cudaGetErrorString(result));
assert(result == cudaSuccess);
}
#endif
return result;
}
__global__ void div_elementwise(const u32 n, float* source, float* dest, const float value)
{
for (u32 i = GET_OFFSET(); i < n; i += GET_STRIDE())
{
dest[i] = source[i] * value;
}
}
float check_eq(float* dev, float* host, u32 len)
{
float sum = 0.0f;
for (u32 i = 0; i < len; ++i)
{
if (dev[i]!=host[i])
{
//printf("diff %d %f %fn", i, dev[i], host[i]);
//break;
float diff = abs((host[i]-dev[i])/host[i]);
sum += diff;
if (diff > (TOL))
printf("diff %d %fn", i, diff);
}
}
printf("%fn", sum);
return sum;
}
void div_host(float* a, float v, u32 len)
{
for (u32 i = 0; i < len; ++i)
{
a[i]=a[i]*v;
}
}
int main()
{
u32 len = sizeof(array)/sizeof(float);
printf("array len = %dn", len);
for (int i =0; i < len; i++) if (isnan(array[i])) {printf("nan value at %dn",i); return -1;}
float* adev, *adevcublas, *d_zero;
float* ahost = (float*) malloc(len * sizeof(float));
checkCuda(cudaMalloc(&adev, len * sizeof(float)));
checkCuda(cudaMalloc(&adevcublas, len * sizeof(float)));
checkCuda(cudaMalloc(&d_zero, len * sizeof(float)));
memcpy(ahost, &array[0], len * sizeof(float));
checkCuda(cudaMemcpy(adev, ahost, len * sizeof(float), cudaMemcpyHostToDevice));
checkCuda(cudaMemcpy(adevcublas, ahost, len * sizeof(float), cudaMemcpyHostToDevice));
checkCuda(cudaMemset(d_zero, 0, len*sizeof(float)));
float alpha = 1/2494.f;
printf("%fn", alpha);
div_host(ahost, alpha, len);
u32 tb = 256;
div_elementwise<<<((len + tb - 1) / tb),tb>>>(len, adev, adev, alpha);
float* r = (float*) malloc(len * sizeof(float));
checkCuda(cudaMemcpy(r, adev, len * sizeof(float), cudaMemcpyDeviceToHost));
check_eq(r,ahost,len);

cublasHandle_t ch;
cublasCreate(&ch);
float* r0 = (float*) malloc(len * sizeof(float));
cublasStatus_t stat = cublasSaxpy(ch, len, &alpha, adevcublas, 1, d_zero, 1);
if (stat != CUBLAS_STATUS_SUCCESS) {std::cout << "CUBLAS error: " << (int)stat << std::endl; return 1;}
checkCuda(cudaMemcpy(r0, d_zero, len * sizeof(float), cudaMemcpyDeviceToHost));

check_eq(r0,ahost,len);
free(r);
free(r0);
free(ahost);
cudaFree(adev);
return 0;
}
相关文章：