Eigen::VectorXd::operator += 似乎比循环 std::vector 慢 ~69%
Eigen::VectorXd::operator += seems ~69% slower than looping through a std::vector
下面的代码（需要 google 基准测试）填充两个向量并将它们相加，将结果存储在第一个向量中。对于我使用 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 结尾）。
FWIW：我添加了 /Qvec-report:2 并且编译器报告：
[...]\Core\AssignEvaluator.h(415) : info C5002: loop not vectorized due to reason '1305'
和reason 1305表示“类型信息不足”。
我有根据的猜测是，Eigen 努力使用内在函数（这里 _mm256_add_pd）会适得其反，并且会混淆编译器。让编译器做它的事情（自动矢量化）似乎是一个更好的主意。我是否遗漏了什么或者这可以被视为 Eigen 错误（错过优化机会）？
TL;DR: 问题主要来自constant loop bound 而不是直接来自Eigen。实际上，在第一种情况下，Eigen 将向量的大小存储在向量属性中，而在第二种情况下，您明确使用常量 N.
聪明的编译器可以使用此信息更积极地展开循环，因为他们知道 N 相当大。展开带有小 N 的循环不是一个好主意，因为代码会更大并且必须由处理器读取。如果代码尚未加载到 L1 缓存中，则必须从其他缓存、RAM 甚至最坏情况下的存储设备中加载。增加的延迟通常大于执行具有小展开因子的顺序循环。这就是编译器并不总是展开循环的原因（至少展开因子不是很大）。
内联在这段代码中也起着重要的作用。事实上，如果函数是内联的，编译器可以 传播常量 并知道向量的大小，使其能够通过更积极地展开循环来进一步优化代码。但是，如果函数未内联，则编译器无法知道循环边界。聪明的编译器仍然可以生成条件算法来优化小循环和大循环，但这会使程序更大并引入小的开销。当代码可以向量化但循环边界未知或者编译时别名未知时（生成的变体数量很快就会很大，因此代码大小也会很大），像 ICC 和 Clang 这样的编译器确实会生成不同的代码替代方案。
请注意，内联函数可能还不够，因为常量传播可能会被处理运行时定义的变量或非内联函数调用的复杂条件所捕获。或者，持续传播的质量可能不足以满足目标示例的要求。
最后，别名 在编译器生成此代码中的 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 指令显然是一个 次优策略 并且可能是一个 编译器 issue/bug。实际上，生成代码的执行速度肯定会受到像您这样的 Intel 处理器上的存储指令的限制。您的处理器每个周期可以执行 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）在此基准测试中的表现非常不同。
下面的代码（需要 google 基准测试）填充两个向量并将它们相加，将结果存储在第一个向量中。对于我使用 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 结尾）。
FWIW：我添加了 /Qvec-report:2 并且编译器报告：
[...]\Core\AssignEvaluator.h(415) : info C5002: loop not vectorized due to reason '1305'
和reason 1305表示“类型信息不足”。
我有根据的猜测是，Eigen 努力使用内在函数（这里 _mm256_add_pd）会适得其反，并且会混淆编译器。让编译器做它的事情（自动矢量化）似乎是一个更好的主意。我是否遗漏了什么或者这可以被视为 Eigen 错误（错过优化机会）？
TL;DR: 问题主要来自constant loop bound 而不是直接来自Eigen。实际上，在第一种情况下，Eigen 将向量的大小存储在向量属性中，而在第二种情况下，您明确使用常量 N.
聪明的编译器可以使用此信息更积极地展开循环，因为他们知道 N 相当大。展开带有小 N 的循环不是一个好主意，因为代码会更大并且必须由处理器读取。如果代码尚未加载到 L1 缓存中，则必须从其他缓存、RAM 甚至最坏情况下的存储设备中加载。增加的延迟通常大于执行具有小展开因子的顺序循环。这就是编译器并不总是展开循环的原因（至少展开因子不是很大）。
内联在这段代码中也起着重要的作用。事实上，如果函数是内联的，编译器可以 传播常量 并知道向量的大小，使其能够通过更积极地展开循环来进一步优化代码。但是，如果函数未内联，则编译器无法知道循环边界。聪明的编译器仍然可以生成条件算法来优化小循环和大循环，但这会使程序更大并引入小的开销。当代码可以向量化但循环边界未知或者编译时别名未知时（生成的变体数量很快就会很大，因此代码大小也会很大），像 ICC 和 Clang 这样的编译器确实会生成不同的代码替代方案。
请注意，内联函数可能还不够，因为常量传播可能会被处理运行时定义的变量或非内联函数调用的复杂条件所捕获。或者，持续传播的质量可能不足以满足目标示例的要求。
最后，别名 在编译器生成此代码中的 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 指令显然是一个 次优策略 并且可能是一个 编译器 issue/bug。实际上，生成代码的执行速度肯定会受到像您这样的 Intel 处理器上的存储指令的限制。您的处理器每个周期可以执行 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）在此基准测试中的表现非常不同。