59

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

这让我很惊讶，因为大多数雅可比是从 0 和 1 向前传播的，所以工作量应该是函数的 2-4 倍，而不是 10-12 倍。为了模拟大量的 jacobian 工作是什么样的，我做了一个超级最小的例子，只有一个点积（而不是真实情况下的 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 根本不做任何优化。

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

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

4

3 回答 3

74

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

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

如果您想在本地禁用显式矢量化并且不怕弄乱 Eigen 的内部，您可以通过专门针对此类型引入一个DontVectorize选项Matrix并禁用矢量化：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

于 2018-08-31T11:37:50.767 回答
39

我很失望地发现，在没有启用快速数学的情况下，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：那是在 Eigen 上。大概他们对矩阵运算进行了矢量化，这对于编译器来说更难看穿。使用这些选项正确优化了一个简单的循环：https ://godbolt.org/z/OppEhY

于 2018-08-31T11:23:59.463 回答
12

强制编译器优化 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 不需要并且需要一些帮助。

于 2018-08-31T11:00:41.157 回答