4

我正在处理在 CUDA 中准确计算零阶 I0 的修正贝塞尔函数的问题。

很长一段时间以来，我一直根据论文使用有理切比雪夫近似

JM Blair，“修正贝塞尔函数 I_0(x) 和 I_1(x) 的有理切比雪夫近似”，数学。计算机，卷。28，名词。126，第 581-583 页，1974 年 4 月。

与 Matlab 提供的结果相比，它给出了 1e-29 量级的平均误差。不幸的是，对于我正在开发的新应用程序来说，这种看似高的准确性已经不够了。

Matlab 使用 DE Amos 开发的 Fortran 例程

Amos, DE，“复杂参数和非负阶贝塞尔函数的子程序包”，桑迪亚国家实验室报告，SAND85-1018，1985 年 5 月。

Amos，DE，“用于复杂参数和非负阶贝塞尔函数的便携式软件包”，反式。数学。软件，1986 年。

可从netlib/amos网站下载。

有一些方法可以在 C/C++ 代码中使用这些 Fortran 例程，方法是将它们编译到库文件中，然后使用 C/C++ 包装器（例如参见netlib_wrapping）。我想知道是否有任何方法可以从这些 Fortran 例程中生成设备功能，然后由 CUDA 内核调用）。

有关问题的更多详细信息

我有两个代码，一个是用 Matlab 编写的，一个是用 CUDA 编写的。两者都通过三个步骤进行操作：

1)通过修改后的贝塞尔函数 I0 和数据的零填充进行缩放；

2）快速傅里叶变换；

3）插值。

我将两者与“精确”结果进行比较：作为步骤 3) 的输出，Matlab 给出的相对均方根误差为 1e-10%，而 CUDA 为 1e-2%，所以我开始调查原因。

两个代码的第一步之间的均方根差，即100*sqrt(sum(abs(U_Matlab_step_1-U_CUDA_step_1).^2))/sqrt(sum(abs(U_Matlab_step_1).^2))，是0%（mean(mean(abs(U_Matlab-U_CUDA)))=6e-29而不是）所以我会说它很好。不幸的是，当我转到第 2 步时，错误上升到2e-4%. 最后，如果我将 CUDA 的步骤 2) 与 Matlab 的步骤 1) 的输出一起输入，那么步骤 2) 的 rms 误差变为1e-14%，这让我认为不准确的来源是由于第一步，即计算修正的贝塞尔函数。

对于本次讨论的有趣发展

查看NVIDIA 开发者专区论坛

4

2 回答 2

4

我想知道这是否可以归因于浮点运算之间的精度差异。

有几件事要检查

  1. Cuda 5 添加了一些新的三角函数，它们可能更好地匹配您的计算格式。此外，我认为自第 4 版以来的 CUDA 数学库具有一些贝塞尔函数，尽管我不确定这是否属实或它们与您的问题有多大关系。
  2. 能写个串口CPU版本来测试吗？这将告诉您您的精度问题是否是由于优化导致的，例如使用 64 位与 80 位表示的数字。关闭优化后，您的计算机将主要处理 80 位表示（也许 matlab 会这样做），而打开数学优化后，您的编译器可能会处理不太准确的 64 位表示。这与 x87 和 SSE 之间的差异有关。
  3. 不同的计算能力硬件的精度略有不同。例如，compute 2.0 执行的 FMA 更准确，更接近优化的 x86。
  4. 是否有物理理由认为 Matlab 是正确的？可能是您的算法在 Matlab 过冲时低于结果。如果 CUDA 对操作进行分组，而 Matlab 没有，则可能会发生这种情况。
  5. 如果必须，必须重新创建 Matlab 结果，您可以尝试通过将输出与不同的舍入技巧匹配来调整代码的每个步骤。见表。

圆桌会议

addition       | x + y        | __dadd_[rn|rz|ru|rd](x, y)
multiplication | x * y        | __dmul_[rn|rz|ru|rd](x, y)
Fused-Mult-Add | fma(x, y, z) | __fma_[rn|rz|ru|rd](x, y, z)
reciprocal     | 1.0 / x      | __drcp_[rn|rz|ru|rd](x)
division       | x / y        | __ddiv_[rn|rz|ru|rd](x, y)
square root    | sqrt(x)      | __dsqrt_[rn|rz|ru|rd](x)

mode | interpretation
rn   | round to nearest, ties to even
rz   | round towards zero
ru   | round towards +∞
rd   | round towards -∞

来自http://developer.download.nvidia.com/assets/cuda/files/NVIDIA-CUDA-Floating-Point.pdf

于 2013-01-02T21:30:32.527 回答
0

我找到了一个介绍性的技术演讲来回答你的问题。这是PDF的链接。所以是的，这是可能的，但是我无法使用上述脚本将旧版 fortran 代码转换为 CUDA C，也许可以直接联系开发人员。

于 2013-01-02T15:19:03.523 回答