替换极慢的pow()函数

Question

我们有一个CFD求解器，运行模拟时发现在某些机器上运行特别慢，但在其他机器上没有这个问题。使用Intel VTune，发现以下代码行是问题所在（Fortran语言）：

RHOV= RHO_INF*((1.0_wp - COEFF*EXP(F0)))**(1.0_wp/(GAMM - 1.0_wp))

通过使用 VTune 调试工具，问题被追踪到 call pow 汇编行，并且在跟踪堆栈时，显示它正在使用 __slowpow()。经过一些搜索，发现了这个页面 抱怨同样的事情。

在 libc 版本为 2.12 的计算机上，模拟需要 18 秒。在 libc 版本为 2.14 的计算机上，模拟只需 0 秒。

基于前面提到的网页信息，当底数接近 1.0 时，会出现问题。因此，我们进行了另一个简单的测试，在 pow() 之前将底数乘以任意数字，然后在调用 pow() 后将结果除以数字的指数幂。这将运行时间从 18 秒降至 0 秒，即使在 libc 2.12 上也是如此。

但是，在所有使用 a**b 的代码中都这样做是不切实际的。有什么方法可以替换 libc 中的 pow() 函数吗？例如，我希望 Fortran 编译器生成的汇编行 call pow 调用我们编写的自定义 pow() 函数，该函数执行缩放、调用 libc 的 pow() 函数，然后除以缩放。如何创建一个中间层对编译器透明呢？ 编辑

为澄清起见，我们正在寻找类似于以下的东西（伪代码）：

double pow(a,b) {
   a *= 5.0
   tmp = pow_from_libc(a,b)
   return tmp/pow_from_libc(5.0, b)
}

是否可以从libc中加载pow并在我们的自定义函数中重命名它以避免命名冲突？如果customPow.o文件可以将libc中的pow重命名，那么如果libc仍然需要用于其他事情，会发生什么情况？这会在customPow.o和libc中导致pow的命名冲突吗？

- tpg2114

---

好老的Fortran！不过问题很有趣 +1 - Austin Henley

Answer 1

好的，请稍等。图书馆不是为了愚弄你而调用__slowpow()函数；它是因为相信对于你所输入的值（在这种情况下，底数接近1，指数数量级为1），需要额外的精度才能给出准确的结果。如果您关心此计算的准确性，在尝试解决问题之前，应该理解其中的原因并确定其是否重要。例如，对于（比如）非常大的负F0值，整个计算过程都可以安全地舍入为1；或者可能不行，这取决于以后对该值的处理方式。如果您需要将1.d0减去这个结果，那么您会想要那些额外的精度。

- Jonathan Dursi

---

这是确实的。但是，至少在我们的情况下，我们的代码是有维度的，所以只有在计算可视化或某种后处理时，我们才会有一个接近于1的基数，因此，接近正确答案的1e-15并不是非常重要。我进行了比较，看看我们失去了多少，误差约为1e-13，对于我们的二阶精度代码来说，比离散化误差小，因此用稍微不太准确的pow()替换所有pow()对我们来说是安全的。 - tpg2114

Answer 2

只需编写您自己的pow函数，将.o文件放入一个静态库归档文件libmypow.a中，存放在链接器的库路径中，并在链接时传递-lmypow参数。

- R.. GitHub STOP HELPING ICE

---

这样做是否允许自定义的 pow 函数调用 libc 中的 pow 呢？如果需要，此自定义 pow 函数将只是对基数进行缩放，然后调用 libc 中的 pow 函数，最后再进行还原。这似乎会导致一些命名冲突。 - tpg2114

如果您使用动态链接，可以使用“dlsym”技巧来实现所需的行为，但这种方法不够稳定。更好的方法是，如果您只需要在使用GNU链接器的系统上工作，则可以使用ld的“--wrap”选项（gcc可以通过-Wl,--wrap,pow将其传递给ld）。然后，您将__wrap_pow放入libmypow.a中，并在需要使用libc pow的地方调用__real_pow，所有事情都应该很好。 - R.. GitHub STOP HELPING ICE

Answer 3

pow(a,b)相当于exp(b*ln(a))，也许这个替换对你有用。

- Mark Ransom

---

这可能会规避调用缓慢的问题，但我们正在寻找一种方法，在不改变现有代码库的情况下，实质上替换Fortran中**运算符生成的函数调用。如果可能的话。 - tpg2114

那么请链接您自己的pow()版本，使用这个等式。 - Hans Passant

这里提供了针对1.0000000000000020^1.5的不同结果：使用call pow得到1.0000000000000031，使用-ffast-math得到1.0000000000000029，而使用exp(b*ln(a))则得到1.5000000000000013。 - steabert

Answer 4

我亲自测试了一下，确实如果我从你提供的链接页面编译测试程序，它会在汇编代码中使用call pow。然而，如果使用优化选项-ffast-math进行编译，就不会调用pow函数，但结果略有不同。

- steabert