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

Posted 2023-02-21

【中文标题】使用特定输入的 cuda/cublas 简单内核中的数值错误

【英文标题】：Numerical error in cuda/cublas simple kernel using particular input

【发布时间】：2015-04-16 14:45:31

【问题描述】：

我正在使用 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 here（您可以使用 nvcc mcve.cu -lcublas 编译它。正如您所看到的两种情况下的结果都是完全错误的：主机结果完全不同，我试图弄清楚发生了什么我在代码中没有看到任何错误，但也许我应该尝试使用双精度而不是浮点数，看看会发生什么。 对这种情况有何看法？提前致谢！

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

【参考方案1】：

有趣的 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。一般来说，我不希望这些结果相同。

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

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

if 测试不总是正确的吗？也许您的意思是1e-12，尽管这仍然存在缺陷。在浮点比较上寻找一个固定的错误阈值应该缩放到被比较的数字的大小。 float 数量仅包含大约 6-7 个准确的十进制数字。 （而且对这些错误求和也很麻烦。）

这是一个完整的代码，修复了上述问题，并为所有比较（主机内核和主机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: %s\n", 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 %f\n", i, dev[i], host[i]);
                        //break;
                        float diff = abs((host[i]-dev[i])/host[i]);
                        sum += diff;
                        if (diff > (TOL))
                                printf("diff %d %f\n", i, diff);
                
        
        printf("%f\n", 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 = %d\n", len);
        for (int i =0; i < len; i++) if (isnan(array[i])) printf("nan value at %d\n",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("%f\n", 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;