为什么 C++ 编译器不做更好的常量折叠？

Question

我正在研究加速大部分 C++ 代码的方法，这些代码具有用于计算 jacobi 的自动导数。这涉及在实际残差中做一些工作，但大部分工作（基于分析的执行时间）是在计算雅可比矩阵。

这让我很惊讶，因为大多数雅可比矩阵都是从 0 和 1 向前传播的，所以工作量应该是函数的 2-4 倍，而不是 10-12 倍。为了模拟大量的 jacobian 工作是什么样的，我做了一个超级小的例子，只有一个点积（而不是真实情况下的 sin、cos、sqrt 等），编译器应该能够优化为单个 return 值：

#include <Eigen/Core>
#include <Eigen/Geometry>

using Array12d = Eigen::Matrix<double,12,1>;

double testReturnFirstDot(const Array12d& b)
{
    Array12d a;
    a.array() = 0.;
    a(0) = 1.;
    return a.dot(b);
}

应该和

一样

double testReturnFirst(const Array12d& b)
{
    return b(0);
}

我失望地发现，如果没有启用快速数学运算，GCC 8.2、Clang 6 或 MSVC 19 都无法对矩阵全为 0 的朴素点积进行任何优化。即使使用快速数学 (https://godbolt.org/z/GvPXFy)，GCC 和 Clang 中的优化也很差（仍然涉及乘法和加法），MSVC 根本不做任何优化。

我没有编译器方面的背景，但这有什么原因吗？我相当确定，在很大一部分科学计算中，能够做得更好的常量 propagation/folding 会使更多的优化变得明显，即使常量折叠本身并没有导致加速。

虽然我对解释为什么不在编译器端完成这件事很感兴趣，但我也对在面对这些类型的模式时我可以在实际方面做些什么来使我自己的代码更快感兴趣.

Answer 1

强制编译器优化 0 和 1 的乘法的一种方法是手动展开循环。为简单起见，我们使用

#include <array>
#include <cstddef>
constexpr std::size_t n = 12;
using Array = std::array<double, n>;

然后我们可以使用折叠表达式（或递归，如果它们不可用）实现一个简单的 dot 函数：

<utility>
template<std::size_t... is>
double dot(const Array& x, const Array& y, std::index_sequence<is...>)
{
    return ((x[is] * y[is]) + ...);
}

double dot(const Array& x, const Array& y)
{
    return dot(x, y, std::make_index_sequence<n>{});
}

现在让我们来看看你的函数

double test(const Array& b)
{
    const Array a{1};    // = {1, 0, ...}
    return dot(a, b);
}

使用 -ffast-math gcc 8.2 produces:

test(std::array<double, 12ul> const&):
  movsd xmm0, QWORD PTR [rdi]
  ret

clang 6.0.0 遵循相同的思路：

test(std::array<double, 12ul> const&): # @test(std::array<double, 12ul> const&)
  movsd xmm0, qword ptr [rdi] # xmm0 = mem[0],zero
  ret

例如，对于

double test(const Array& b)
{
    const Array a{1, 1};    // = {1, 1, 0...}
    return dot(a, b);
}

我们得到

test(std::array<double, 12ul> const&):
  movsd xmm0, QWORD PTR [rdi]
  addsd xmm0, QWORD PTR [rdi+8]
  ret

添加。 Clang 展开了一个 for (std::size_t i = 0; i < n; ++i) ... 循环，没有所有这些折叠表达式技巧，gcc 没有并且需要一些帮助。

Answer 2

I was disappointed to find that, without fast-math enabled, neither GCC 8.2, Clang 6 or MSVC 19 were able to make any optimizations at all over the naive dot-product with a matrix full of 0s.

不幸的是，他们别无选择。由于 IEEE 浮点数有符号零，添加 0.0 不是恒等运算：

-0.0 + 0.0 = 0.0 // Not -0.0!

同样，乘以零并不总是得到零：

0.0 * Infinity = NaN // Not 0.0!

因此，编译器根本无法在保持 IEEE 浮点合规性的同时在点积中执行这些常量折叠 - 据他们所知，您的输入可能包含带符号的零 and/or 无穷大。

您将不得不使用 -ffast-math 来获得这些折叠，但这可能会产生不良后果。您可以使用特定标志（来自 http://gcc.gnu.org/wiki/FloatingPointMath）获得更细粒度的控制。根据上面的解释，加上下面两个flag应该可以实现常量折叠：
-ffinite-math-only、-fno-signed-zeros

确实，您可以通过这种方式获得与 -ffast-math 相同的程序集：https://godbolt.org/z/vGULLA。您只放弃带符号的零（可能无关紧要）、NaN 和无穷大。据推测，如果您仍然在代码中生成它们，您会得到未定义的行为，因此请权衡您的选择。

至于为什么你的例子即使使用 -ffast-math 也没有得到更好的优化：那是在 Eigen 上。据推测，他们对矩阵运算进行了矢量化，编译器更难看穿。使用这些选项可以正确优化一个简单的循环：https://godbolt.org/z/OppEhY

Answer 3

这是因为 Eigen 明确地将您的代码矢量化为 3 个 vmulpd、2 个 vaddpd 和剩余 4 个组件寄存器中的 1 个水平缩减（这假设是 AVX，只有 SSE，您将获得 6 个 mulpd 和 5 个 addpd）。 -ffast-math 允许 GCC 和 clang 删除最后 2 个 vmulpd 和 vaddpd（这就是他们所做的），但它们不能真正替换 Eigen 明确生成的剩余 vmulpd 和水平缩减。

那么，如果您通过定义 EIGEN_DONT_VECTORIZE 来禁用 Eigen 的显式矢量化呢？然后你得到了你所期望的 (https://godbolt.org/z/UQsoeH) 但其他代码片段可能会变得更慢。

如果你想在本地禁用显式矢量化并且不怕弄乱 Eigen 的内部，你可以向 Matrix 引入一个 DontVectorize 选项并通过专门化 traits<> 来禁用矢量化这个 Matrix 类型：

static const int DontVectorize = 0x80000000;

namespace Eigen {
namespace internal {

template<typename _Scalar, int _Rows, int _Cols, int _MaxRows, int _MaxCols>
struct traits<Matrix<_Scalar, _Rows, _Cols, DontVectorize, _MaxRows, _MaxCols> >
: traits<Matrix<_Scalar, _Rows, _Cols> >
{
  typedef traits<Matrix<_Scalar, _Rows, _Cols> > Base;
  enum {
    EvaluatorFlags = Base::EvaluatorFlags & ~PacketAccessBit
  };
};

}
}

using ArrayS12d = Eigen::Matrix<double,12,1,DontVectorize>;

完整示例：https://godbolt.org/z/bOEyzv

为什么 C++ 编译器不做更好的常量折叠？

Why don't C++ compilers do better constant folding?

c++

compiler-construction

automatic-differentiation

eigen

ceres-solver