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、/D "NDEBUG"、/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 个 double。但是，对于 std::vector，编译器对循环进行了展开，每次迭代执行 2 次 vaddpd，但是对于 Eigen::VectorXd，它没有这样做。另一个可能重要的区别是，std::vector 的循环对齐到 32 字节（地址以 0x40 = 64 = 2*32 结尾）。

顺便提一下：我已经添加了 /Qvec-report:2，编译器报告如下：

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

而 1305 原因代码 意味着“缺少类型信息”。

我的猜测是，Eigen 使用内部函数（例如 _mm256_add_pd）实现向量化可能会适得其反，混淆编译器。让编译器自行进行自动向量化似乎是一个更好的想法。我是否有所遗漏？或者这可以被视为 Eigen 的缺陷（错失的优化机会）？

- Cassio Neri

Answer 1

TL;DR: 问题主要来自于循环边界常量，而不是直接来自于Eigen。事实上，在第一种情况下，Eigen将向量大小存储在向量属性中，而在第二种情况下，您明确使用常量N。

聪明的编译器可以利用此信息更积极地展开循环，因为它们知道N相当大。展开具有小N的循环是一个坏主意，因为代码会变得更大，必须由处理器读取。如果代码尚未加载到L1缓存中，则必须从其他缓存、RAM甚至最坏的情况下从存储设备加载。增加的延迟通常比执行具有小展开因子的顺序循环更大。这就是为什么编译器不总是展开循环（至少不使用大的展开因子）的原因。 内联也在这段代码中发挥着重要作用。确实，如果函数被内联，编译器可以传播常量并知道向量的大小，使其能够通过更积极地展开循环进一步优化代码。但是，如果函数没有被内联，那么编译器无法知道循环边界。聪明的编译器仍然可以生成条件算法以优化小循环和大循环，但这会使程序变得更大并引入一些开销。像ICC和Clang这样的编译器在代码可以矢量化但循环边界未知或别名也未知时生成不同的代码替代方案（生成的变体数量可能很快就会非常庞大，因此代码大小也会增加）。

请注意，内联函数可能还不够，因为复杂的条件语句涉及运行时定义的变量或非内联函数调用可能会阻碍常量传播。或者，常量传播的质量可能不足以针对目标示例进行优化。

最后，别名在编译器生成SIMD指令（并可能更好地展开循环）的能力中也起着重要作用。实际上，别名经常会阻止使用SIMD指令，而编译器并不总是容易检查别名并相应地生成快速实现。

测试假设

如果基于向量的实现使用存储在向量对象中的循环边界，则MSVC生成的代码在基准测试中未被向量化：尽管函数已内联，但常量未被正确传播。结果代码应该更慢。这是生成的代码：

$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生成的代码在基准测试中将会被向量化并正确展开：编译时的常量帮助编译器生成一个循环展开两次的循环。它通过混合SSE和AVX指令来完成，这是非常令人惊讶的（该点将在下面讨论）。生成的代码应该比原始的Eigen实现快得多。然而，由于意外使用了SSE指令，它可能不如初始的向量实现快。这里是生成的代码：

$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指令，此代码仍然不完美。

或者，可以尝试使用固定大小的Eigen向量以获得更好的性能。

最后，请注意，其他编译器（例如Clang、ICC和GCC）在这个基准测试中的行为非常不同。

- Jérôme Richard

---

非常准确！如果我将add_vectors函数中针对std::vector重载的循环边界改为v2.size()（而不是N），则性能会明显改善，使Eigen库成为最佳选择（197ns x 739ns）。 - Cassio Neri