为什么计算科学家需要实现他们自己的 std::complex 版本?

计算科学 线性代数 C++ 编程范式 复数
2021-12-20 23:01:06

计算科学中许多著名的 C++ 库,例如EigenTrilinosdeal.II使用标准 C++ 模板头库对象std::complex<>, 来表示复杂的浮点数。

在 Jack Poulson 对有关默认构造函数的问题的回答中,他指出“出于多种原因” ,他std::complexElemental中有自己的实现。这些原因是什么?这种方法的优点和缺点是什么?

3个回答

我相信这个讨论已经多次出现在 PETSc 列表中。我的主要原因是:

  1. C++ 标准规定 std::complex 仅为 float、double 和 long double 数据类型定义。因此它不能用于其他数据类型,例如四精度。

  2. 该标准不保证复杂算术的稳定性。

  3. 该标准不保证 std::complex 中的数据存储为实部后跟虚部。这对于与外部库(例如 BLAS 和 LAPACK)的接口至关重要。所有主要实现都是如此,但我希望能够确保它。

  4. 我更喜欢能够直接操纵实部和虚部。std::complex 使这变得不必要地困难。

  5. 我希望最终有一个更通用的版本,它只需要数据类型是一个环而不需要一个字段。这将包括高斯整数。

std::complex<>在我的程序中使用,并且必须为每个新的编译器或编译器升级与编译器标志和解决方法作斗争。我将尝试按时间顺序讲述这些战斗:

  1. 性能测量表明,仅涉及计算复数域的绝对值平方的步骤比之前的 gcc-4.x 的 FFT 花费更多时间。深入研究生成的汇编代码,发现std::norm( ) 以一种避免溢出的方式计算了绝对值 ( ),然后对结果求平方。这个问题可以通过编译标志来解决|z|2|z|-ffast-math
  2. linux(或链接器)上的 intel icc 编译器std::arg在某些配置下编译为非 opt(与特定 gcc 版本的链接兼容性)。这个问题经常出现,所以std::arg不得不替换为atan2(imag(),real()). 但是在编写新代码时很容易忘记这一点。
  3. 该类型std::complex使用不同于内置 C99 复杂类型的调用约定 (=ABI),以及用于较新 gcc 版本的内置 Fortran 复杂类型。
  4. -ffast-math编译标志以意想不到的方式与浮点异常的处理交互发生的情况是编译器将除法拉出循环,从而division by zero在运行时导致异常。这些异常永远不会在循环内部发生,因为由于周围的逻辑,相应的除法没有发生。那个真的很糟糕,因为它是一个与使用浮点异常处理(使用不同的编译标志)的程序分开编译的库,并遇到了这些问题(相应的团队坐在世界的另一端,所以这个问题确实造成了很大的麻烦)。这是通过更加小心地手动进行编译器使用的优化来解决的。
  5. 该库成为程序的一部分,不再使用-ffast-math编译标志。升级到较新的 gcc 版本后,性能下降了一个很大的因素。我还没有详细调查过这个问题,但我担心它与C99 Annex G有关。我不得不承认,我完全被这个奇怪的复数乘法定义弄糊涂了,甚至似乎存在不同的版本,声称其他版本被误导了。我希望-fcx-limited-range编译标志能解决这个问题,因为这个较新的 gcc 版本似乎还有另一个问题-ffast-math
  6. 编译标志使新版本的 gcc 的行为完全不可预测(甚至受到-ffast-math影响)。唯一的解决方法似乎是避免在程序中出现任何.NaNisnanNaNNaN

现在你可能会问我是否打算放弃内置的复杂类型以及std::complex出于这些原因。只要我继续使用 C++,我就会继续使用内置类型。如果 C++ 变得完全无法用于科学计算,我宁愿考虑切换到一种更关注与科学计算相关的问题的语言。

另一个可能的原因可能是历史原因:您自己的版本早于当前公认的标准。

虽然这不是关于复数,但这个错误报告的讨论是关于容器类型的使用。有问题的软件,即 OpenFOAM,早在 C++ 的标准模板库 (STL) 出现之前就已经开发出来了。因此,该软件包含许多数据类型,现在这些数据类型只是从 STL 中获取。

虽然沉没成本谬误通常被认为是一种谬误,但仍然存在沉没成本,这导致许多开发人员维护自己的“标准”类型实现,因为重构他们的代码会占用大量资源。