Eigen::VectorXd::operator += 似乎比通过 std::vector 循环慢约 69%

Question

下面的代码（需要google benchmark）填充两个向量并将它们相加，将结果存储在第一个向量中。对于我使用Eigen::VectorXd 和std::vector 进行性能比较的向量类型：

#include <Eigen/Core>
#include <benchmark/benchmark.h>
#include <vector>

auto constexpr N = 1024u;

template <typename TVector>
TVector generate(unsigned min) {
    TVector v(N);
    for (unsigned i = 0; i < N; ++i)
        v[i] = static_cast<double>(min + i);
    return v;
}

auto ev1 = generate<Eigen::VectorXd>(0);
auto ev2 = generate<Eigen::VectorXd>(N);
auto sv1 = generate<std::vector<double>>(0);
auto sv2 = generate<std::vector<double>>(N);

void add_vectors(Eigen::VectorXd& v1, Eigen::VectorXd const& v2) {
    v1 += v2;
}

void add_vectors(std::vector<double>& v1, std::vector<double> const& v2) {
    for (unsigned i = 0; i < N; ++i)
        v1[i] += v2[i];
}

static void eigen(benchmark::State& state) {
    for (auto _ : state) {
        add_vectors(ev1, ev2);
        benchmark::DoNotOptimize(ev1);
    }
}

static void standard(benchmark::State& state) {
    for (auto _ : state) {
        add_vectors(sv1, sv2);
        benchmark::DoNotOptimize(sv1);
    }
}

BENCHMARK(standard);
BENCHMARK(eigen);

我在 Intel Xeon E-2286M @2.40Ghz 上运行它，使用 Eigen 3.3.9、MSVC 16.11.2 以及这些相关的编译器开关 /GL、/Gy、/O2、@ 987654328@、/Oi 和 /arch:AVX。典型的输出如下所示：

Run on (16 X 2400 MHz CPU s)
CPU Caches:
  L1 Data 32K (x8)
  L1 Instruction 32K (x8)
  L2 Unified 262K (x8)
  L3 Unified 16777K (x1)
--------------------------------------------------
Benchmark           Time           CPU Iterations
--------------------------------------------------
standard           99 ns        100 ns    7466667
eigen             169 ns        169 ns    4072727

这似乎表明在 std::vector 上的操作比在 Eigen::VectorXd 上快约 69%。在反汇编中，紧密的循环如下所示：

// For Eigen::VectorXd
00007FF672221A11  vmovupd     ymm0,ymmword ptr [rcx+rax*8]  
00007FF672221A16  vaddpd      ymm1,ymm0,ymmword ptr [r8+rax*8]  
00007FF672221A1C  vmovupd     ymmword ptr [r8+rax*8],ymm1  
00007FF672221A22  add         rax,4  
00007FF672221A26  cmp         rax,rdx  
00007FF672221A29  jge         eigen+0C7h (07FF672221A37h)  
00007FF672221A2B  mov         rcx,qword ptr [rsp+48h]  
00007FF672221A30  mov         r8,qword ptr [rsp+58h]  
00007FF672221A35  jmp         eigen+0A1h (07FF672221A11h)  

// For std::vector
00007FF672221B40  vmovups     ymm1,ymmword ptr [rax+rdx-20h]  
00007FF672221B46  vaddpd      ymm1,ymm1,ymmword ptr [rax+rcx-20h]  
00007FF672221B4C  vmovups     ymmword ptr [rax+rcx-20h],ymm1  
00007FF672221B52  vmovups     ymm1,ymmword ptr [rax+rdx]  
00007FF672221B57  vaddpd      ymm1,ymm1,ymmword ptr [rax+rcx]  
00007FF672221B5C  vmovups     ymmword ptr [rax+rcx],ymm1  
00007FF672221B61  lea         rax,[rax+40h]  
00007FF672221B65  sub         r8,1  
00007FF672221B69  jne         standard+0C0h (07FF672221B40h)  

可以注意到，两者都使用vaddpd 来一次添加4 个doubles。但是，对于std::vector，编译器展开循环以在每次迭代中执行 2 个vaddpd，但对于Eigen::VectorXd 却没有这样做。另一个潜在的重要区别是std::vector 的循环对齐到 32 个字节（地址以 0x40 = 64 = 2*32 结尾）。

FWIW：我添加了 /Qvec-report:2 并且编译器报告：

[...]\Core\AssignEvaluator.h(415) : info C5002: loop not vectorized due to reason '1305'

而reason 1305 表示“类型信息不足”。

我有根据的猜测是，Eigen 使用内在函数（此处为 _mm256_add_pd）的努力会适得其反，并使编译器感到困惑。让编译器做它的事情（自动矢量化）似乎是一个更好的主意。我是否遗漏了什么，或者这可以被认为是一个 Eigen 错误（错过了优化机会）？

Answer 1

TL;DR: 问题主要来自constant loop bound，而不是直接来自Eigen。实际上，在第一种情况下，Eigen 将向量的大小存储在向量属性中，而在第二种情况下，您显式使用常量 N。

聪明的编译器可以使用这些信息更积极地展开循环，因为他们知道N 相当大。使用小的N 展开循环是个坏主意，因为代码会更大并且必须由处理器读取。如果代码尚未加载到 L1 缓存中，则必须从其他缓存、RAM 甚至在最坏情况下的存储设备中加载。增加的延迟通常大于执行具有较小展开因子的顺序循环。这就是为什么编译器并不总是展开循环的原因（至少没有很大的展开因子）。

内联在这段代码中也扮演着重要的角色。事实上，如果函数是内联的，编译器可以传播常量并知道向量的大小，使其能够通过更积极地展开循环来进一步优化代码。但是，如果函数没有内联，那么编译器就无法知道循环边界。聪明的编译器仍然可以生成条件算法来优化小循环和大循环，但这会使程序更大并引入少量开销。当代码可以向量化但循环边界未知时，或者在编译时不知道别名时（生成的变体的数量可能很快很大，因此代码大小），ICC 和 Clang 等编译器确实会生成不同的代码替代方案。

请注意，内联函数可能还不够，因为常量传播可能会被处理运行时定义的变量或非内联函数调用的复杂条件捕获。或者，常数传播的质量可能不足以满足目标示例。

最后，aliasing 在编译器在此代码中生成 SIMD 指令（并可能更好地展开循环）的能力方面也发挥着关键作用。事实上，别名通常会阻止 SIMD 指令的使用，而且编译器检查别名并相应地生成快速实现并不总是那么容易。

检验假设

如果基于向量的实现使用存储在向量对象中的循环绑定，则 MSVC 生成的代码在基准测试中未向量化：尽管内联了功能。生成的代码应该慢得多。这是generated code：

$LL24@standard:
        vmovsd  xmm0, QWORD PTR [r9+rcx*8]
        vaddsd  xmm1, xmm0, QWORD PTR [r8+rcx*8]
        vmovsd  QWORD PTR [r8+rcx*8], xmm1
        mov     rax, QWORD PTR std::vector<double,std::allocator<double> > sv1+8
        inc     edx
        sub     rax, QWORD PTR std::vector<double,std::allocator<double> > sv1
        sar     rax, 3
        mov     ecx, edx
        cmp     rcx, rax
        jb      SHORT $LL24@standard

如果基于 Eigen 的实现使用常量循环绑定，则 MSVC 生成的代码在基准测试中矢量化和正确展开：编译时常量有助于编译器生成循环展开 2 次。它通过混合 SSE 和 AVX 指令来做到这一点，这非常令人惊讶（这一点将在下面讨论）。生成的代码应该比原始的 Eigen 实现快得多。但是，由于意外使用 SSE 指令，它可能不如初始向量实现快。这是generated code：

$LL24@eigen:
        vmovupd xmm1, XMMWORD PTR [rdx+rcx-16]
        vaddpd  xmm1, xmm1, XMMWORD PTR [rcx-16]
        vmovupd xmm2, XMMWORD PTR [rcx+rdx]
        vmovupd XMMWORD PTR [rcx-16], xmm1
        vaddpd  xmm1, xmm2, XMMWORD PTR [rcx]
        vmovupd XMMWORD PTR [rcx], xmm1
        vmovups ymm1, YMMWORD PTR [rdx+rcx+16]
        vaddpd  ymm1, ymm1, YMMWORD PTR [rcx+16]
        vmovups YMMWORD PTR [rcx+16], ymm1
        lea     rcx, QWORD PTR [rcx+64]
        sub     rax, 1
        jne     SHORT $LL24@eigen

补充说明

值得注意的是，为非内联版本生成的代码使用了非常低效的标量代码（通常是由于 N 未知并且指针别名可能是可能的）。

在您的情况下，在这样的循环中混合 SSE 和 AVX 指令显然是一个次优策略，并且可能是一个编译器问题/错误。事实上，生成的代码的执行速度肯定会受到像您这样的英特尔处理器上的存储指令的限制。您的处理器每个周期可以执行 1 个存储指令，每个周期可以执行 2 个加载指令，并且每个周期可以计算 2 个向量化加法。它每个周期最多可以执行 6 条微指令（来自 5 条独立指令，可能还有 4 条缓存的附加指令）。因此，生成的混合 SSE 和 AVX 的代码每次迭代至少需要 3 个周期。同时，最初的基于向量的实现可以在 2 个周期内执行 4 个加载、2 个存储、2 个加法和 3 个其他指令，如 lea/sub/branch（由于实际微指令端口等复杂的硬件，实际上可能执行 3 个）调度，微指令缓存）。但是，请注意，编译器参数并未指定针对您的特定处理器架构（即 Intel Coffee Lake）优化代码。尽管如此，我仍然非常怀疑混合 SSE 和 AVX 代码是否会显着提升 AMD 处理器（或任何主流 x86 处理器）的性能。或者，我可能是因为 MSVC 在这种情况下无法完全检测到没有别名。

要消除妨碍代码矢量化和循环展开的大多数别名问题，可以使用OpenMP SIMD 指令（例如#pragma omp simd）。 MSVC 使用标志/openmp:experimental 实验性地支持这一点。这是生成的代码：

void add_vectors(Eigen::VectorXd& v1, Eigen::VectorXd const& v2) {
    #pragma omp simd
    for (unsigned i = 0; i < N; ++i)
        v1[i] += v2[i];
}

MSVC 出人意料地生成了一个只有 SSE 指令的汇编代码，但是如果你启用了 AVX2，那么它会生成一个相对好的代码：

$LL26@eigen:
        mov     rcx, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev1
        lea     rdx, QWORD PTR [rdx+128]
        mov     rax, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev2
        vmovupd ymm0, YMMWORD PTR [rdx+rcx-192]
        vaddpd  ymm0, ymm0, YMMWORD PTR [rdx+rax-192]
        vmovupd YMMWORD PTR [rdx+rcx-192], ymm0
        mov     rcx, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev1
        mov     rax, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev2
        vmovupd ymm0, YMMWORD PTR [rdx+rcx-160]
        vaddpd  ymm0, ymm0, YMMWORD PTR [rdx+rax-160]
        vmovupd YMMWORD PTR [rdx+rcx-160], ymm0
        mov     rcx, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev1
        mov     rax, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev2
        vmovupd ymm0, YMMWORD PTR [rdx+rcx-128]
        vaddpd  ymm0, ymm0, YMMWORD PTR [rdx+rax-128]
        vmovupd YMMWORD PTR [rdx+rcx-128], ymm0
        mov     rcx, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev1
        mov     rax, QWORD PTR Eigen::Matrix<double,-1,1,0,-1,1> ev2
        vmovupd ymm0, YMMWORD PTR [rdx+rcx-96]
        vaddpd  ymm0, ymm0, YMMWORD PTR [rdx+rax-96]
        vmovupd YMMWORD PTR [rdx+rcx-96], ymm0
        sub     r8, 1
        jne     $LL26@eigen

由于意外的无用 mov 指令，此代码仍然不完美。

或者，也可以使用fixed-size Eigen vectors 以获得更好的性能。

最后，请注意其他编译器（如 Clang、ICC 和 GCC）在此基准测试中的表现非常不同。