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

Question

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

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

#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 根本不进行任何优化。

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

虽然我对解释为什么没有在编译器端完成此操作感兴趣，但我也对在面对这些模式时我可以在实际方面做些什么来使我自己的代码更快感兴趣。

Answer 1

这是因为 Eigen 在剩余的 4 个组件寄存器中将您的代码显式矢量化为 3 个 vmulpd、2 个 vaddpd 和 1 个水平缩减(这假设是 AVX，只有 SSE 你会得到 6 个 mulpd 和 5 个 addpd(。有了-ffast-mathGCC和clang可以删除最后2个vmulpd和vaddpd(这就是他们所做的(，但它们不能真正取代Eigen明确生成的剩余vmulpd和水平约简。

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

如果你想在本地禁用显式矢量化，并且不怕弄乱 Eigen 的内部，你可以引入一个DontVectorize选项来Matrix和禁用矢量化，方法是专门针对此Matrix类型进行traits<>：

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

Answer 2

我很失望地发现，如果没有启用快速数学，GCC 8.2、Clang 6 或 MSVC 19 都无法在充满 0 的矩阵的朴素点积上进行任何优化。

不幸的是，他们别无选择。由于 IEEE 浮点数具有带符号的零，因此添加0.0不是标识操作：

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

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

0.0 * Infinity = NaN // Not 0.0!

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

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

事实上，你会得到与-ffast-math相同的程序集：https://godbolt.org/z/vGULLA。你只放弃有符号的零(可能无关紧要(、NaN 和无穷大。据推测，如果你仍然在代码中生成它们，你会得到未定义的行为，所以权衡你的选择。

---

至于为什么即使使用-ffast-math您的示例也没有更好地优化：那是在本征上。据推测，它们对矩阵操作进行了矢量化处理，这对于编译器来说更难看穿。使用以下选项对简单循环进行适当优化： https://godbolt.org/z/OppEhY

Answer 3

强制编译器优化乘以 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-mathgcc 8.2 产生：

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没有，需要一些帮助。