为什么 CUDA 内核中的 atan2 产生的值与 C/fortran 略有不同？

Question

请帮我解决这个问题：

对于 x = -6.5799015957503127E+02 和 y = -4.6102597302044005E+03

// in C:
atan2(x,y) = -2.99982704867761151845684253203217e+00

// in Fortran:
atan2(x,y) = -2.99982704867761151845684253203217D+00

// But atan2 called in CUDA kernel is:
atan2(x,y) = -2.99982704867761107436763268196955E+00
//                             ^^^^^^^^^^^^^^^^^^^^^

是的，可能是由于舍入错误，但是，为什么在 Fortran 和 C 中结果相同，而在 CUDA 中略有不同？

我在 CUDA 中需要与 Fortran 和 C 中相同数量的 atan2。如何做到这一点？

Answer 1

请注意，您的 x 和 y 值（以十进制表示）不能完全表示为二进制（浮点）数，在 double 的范围内。下面的程序演示了这一点：

$ cat t556.cu
#include <stdio.h>
#include <math.h>

__global__ void mykernel(double x, double y, double *res){

  *res = atan2(x,y);
}

int main(){

  double h_x = -6.5799015957503127E+02;
  double h_y = -4.6102597302044005E+03;
  double *d_res;
  double my_res = atan2(h_x, h_y);

  cudaMalloc(&d_res, sizeof(double));
  mykernel<<<1,1>>>(h_x, h_y, d_res);
  double h_res;
  cudaMemcpy(&h_res, d_res, sizeof(double), cudaMemcpyDeviceToHost);

  printf("x = %.32lf\n", h_x);
  printf("y = %.32lf\n\n", h_y);
  printf("hst = %.32lf\n", my_res);
  printf("dev = %.32lf\n\n", h_res);

  printf("hst bin = 0x%lx\n", *(reinterpret_cast<unsigned long long *>(&my_res)));
  printf("dev bin = 0x%lx\n", *(reinterpret_cast<unsigned long long *>(&h_res)));


  return 0;

}
$ nvcc -arch=sm_20 -o t556 t556.cu
$ ./t556
x = -657.99015957503127083327854052186012
y = -4610.25973020440051186596974730491638

hst = -2.99982704867761151845684253203217
dev = -2.99982704867761107436763268196955

hst bin = 0xc007ffa552ddcff5
dev bin = 0xc007ffa552ddcff4
$

我们看到，当我们根据你在问题中写的内容指定x，然后用大量数字打印出来时，打印输出与代码分配的“大概”值不匹配。

上面的程序还演示了atan2在主机和设备上计算的结果相差1位，在结果尾数的最低位。

参考the CUDA math documentation in the programming guide，我们看到atan2 函数的最大误差为2 ULP（ULP = 排在最后的单位，对于本次讨论，它相当于尾数的最低有效位所代表的单位)。

这意味着atan2 函数（在 CUDA 设备代码中）不能保证产生数字正确的结果（对于整个任意精度），它产生的结果可能与数字正确的结果相差 1 或 2 ULP （自始至终完全任意精度）IEEE-754 实现。这意味着在将atan2 的 CUDA 实现与另一个实现进行比较时，可以合理地假设结果中可能存在 1 或 2 个尾数 LSB 的差异。

如果您要求 CUDA 设备计算的 atan2 结果与另一个实现的 atan2 结果完美匹配（没有尾数位不同），则 CUDA 数学库提供的 atan2 函数将不会对你有用。

那时我能给你的唯一建议是create your own implementation of atan2，在整个过程中使用更基本的浮点运算（可能会在整个过程中选择在 CUDA 数学库中提供 0 ULP 错误的运算，尽管我不是这方面的专家）这将如何详细完成），以及旨在匹配您正在比较的实现的数值方法。

This 也可能是一个信息丰富的阅读。请注意，GNU 实现不一定意味着所有数学函数或什至所有 trig 类型函数的 ULP 错误为 0。例如，请注意 cos 在 IA64 上似乎最多有 1 个 ULP 错误。但是atan2 似乎处于隐含的最低错误级别。

Answer 2

我使用 GPU 架构 sm_20。我项目中的所有函数（sin、cos、sinh、cosh 等）都给出与 C 中相同的结果。那么，为什么 atan2 不呢？？？？ :(

如今，大多数三角函数的实现都精确到略高于 0.5 ULP，这意味着在 99% 的情况下，精确的数学结果只有一个可以返回的可表示浮点近似值（这是最接近实际结果）。

但是，您不应假设任何三角函数都是完美的，即精确到 0.5 ULP，除非您为此属性选择了数学库。这意味着，对于一些数学结果几乎正好在两个可表示的双精度之间的罕见参数，三角函数可能会返回错误的参数（比如距离为 0.507 ULP 的参数，而不是距离为 0.507 ULP 的参数）距离为 0.493 ULP）。

这也意味着两个不同的实现可以返回不同的结果（一个实现可以正确地将结果返回到 0.493 ULP 而另一个错误的结果返回到 0.507 ULP）。

这可能发生在所有三角函数中。你只是碰巧在atan2 遇到了这个问题，但同样的事情也可能发生在sin 或cos 上。在您使用的其中一个库中，atan2 的实现可能不太准确（例如，到 0.52 ULP 而不是 0.505 ULP），这使得问题更容易被注意到。但是除非你在两边都使用正确舍入的库，或者是同一个库（不会被正确舍入，但两边都会产生同样的错误），这种情况会时有发生。 p>

CRlibm 是正确舍入数学库的一个示例，它产生的结果与任何其他正确舍入的数学库相同。 netlib 是一个不错的数学库示例，它经常嵌入到程序中以便它们在任何地方产生相同的结果。