CUDA 浮点精度与 CPU 实现不匹配答案

【问题标题】：CUDA float precision not matching CPU implementationCUDA 浮点精度与 CPU 实现不匹配
【发布时间】：2020-11-10 02:24:34
【问题描述】：

我在 GTX 1080Ti 上使用 CUDA 5.5 计算 3.5 并想计算这个公式：

y = a * a * b / 64 + c * c

假设我有这些参数：

a = 5876
b = 0.4474222958088
c = 664

我同时通过 GPU 和 CPU 进行计算，它们给了我不同的不准确答案：

h_data[0]  = 6.822759375000e+05,
h_ref[0]   = 6.822760000000e+05,
difference = -6.250000000000e-02

h_data 是 CUDA 答案，h_ref 是 CPU 答案。当我将这些插入我的计算器时，GPU 的答案更接近于准确的答案，我怀疑这与浮点精度有关。我现在的问题是，如何让 CUDA 解决方案与 CPU 版本的精度/舍入相匹配？如果我将 a 参数偏移 +/-1，则解决方案匹配，但如果我偏移 c 参数，我仍然会得到 1/16 的差异

这是工作代码：

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

__global__ void test_func(float a, float b, int c, int nz, float * __restrict__ d_out)
{
  float *fdes_out = d_out + blockIdx.x * nz;
  float roffout2 = a * a / 64.f;
  //float tmp = fma(roffout2,vel,index*index);
  for (int tid = threadIdx.x; tid < nz; tid += blockDim.x) {
    fdes_out[tid] = roffout2 * b + c * c;
  }
}

int main (int argc, char **argv) 
{
  // parameters
  float a = 5876.0f, b = 0.4474222958088f;
  int c = 664; 

  int nz = 1;
  float *d_data, *h_data, *h_ref;
  h_data = (float*)malloc(nz*sizeof(float));
  h_ref  = (float*)malloc(nz*sizeof(float));
  
  // CUDA
  cudaMalloc((void**)&d_data, sizeof(float)*nz); 
  dim3 nb(1,1,1); dim3 nt(64,1,1);
  test_func <<<nb,nt>>> (a,b,c,nz,d_data);
  cudaMemcpy(h_data, d_data, sizeof(float)*nz, cudaMemcpyDeviceToHost);
  
  // Reference
  float roffout2 = a * a / 64.f;
  h_ref[0] = roffout2*b + c*c;
  
  // Compare
  printf("h_data[0]  = %1.12e,\nh_ref[0]   = %1.12e,\ndifference = %1.12e\n",
    h_data[0],h_ref[0],h_data[0]-h_ref[0]);

  // Free
  free(h_data); free(h_ref);
  cudaFree(d_data);
  return 0;
}

我只使用-O3 标志进行编译。

【问题讨论】：

另请注意，浮点数的差异仅为 1 LSB。我在 MATLAB 中对此进行了测试，使用浮点数的操作给了我y_float=682276 和双精度数y_double=682275.966221463。请注意，这个“真实”结果的差异小于浮点数的 LSB，这意味着两者几乎相同正确。请注意，GPU 一比 CPU 一更接近双精度答案，加强了@njuffa 的猜测
是的，根据我的测试，当使用-fmad=false编译代码时，差异就消失了。并且 OP 已经声明“GPU 答案更接近确切答案”@njuffa 如果你想提供一个我会投票的答案

标签： cuda precision

【解决方案1】：

一个单精度ulp 的数值差异很小，因为 CUDA 编译器默认应用 FMA 合并，而主机编译器不这样做。可以通过将命令行标志 -fmad=false 添加到 CUDA 编译器驱动程序 nvcc 的调用来关闭 FMA 合并。

FMA 合并是一种编译器优化，其中 FMUL 和从属 FADD 被转换为单个 fused multiply-add 或 FMA 指令。 FMA 指令计算 a*b+c，使得完整的未舍入乘积 a*b 进入与 c 的加法，然后应用最终舍入以产生最终结果。

通常，这具有性能优势，因为执行一条 FMA 指令而不是两条指令 FMUL、FADD，并且所有指令都具有相似的延迟。通常，这也具有准确性优势，因为使用 FMA 消除了一个舍入步骤，并在 a*c 和 c 具有相反符号时防止 subtractive cancellation。

在这种情况下，如 OP 所述，使用 FMA 计算的 GPU 结果比未使用 FMA 计算的主机结果稍微准确一些。使用更高精度的参考，我发现GPU结果的相对误差是-4.21e-8，而主机结果的相对误差是4.95e-8。

【讨论】：