【问题标题】：Is it possible to vectorize myNum += a[b[i]] * c[i]; on x86_64?是否可以矢量化 myNum += a[b[i]] * c[i];在 x86_64 上？
【发布时间】：2011-01-21 22:46:39
【问题描述】：

在 x86_64 上我将使用哪些内在函数来矢量化以下内容（如果甚至可以矢量化）？

double myNum = 0;
for(int i=0;i<n;i++){
    myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double
}

【问题讨论】：

索引在b中的分布是什么？
只是好奇，以下建议是否有助于加快您的代码速度？

标签： x86 x86-64 sse simd vectorization

【解决方案1】：

这是我的尝试，经过全面优化和测试：

#include <emmintrin.h>

__m128d sum = _mm_setzero_pd();
for(int i=0; i<n; i+=2) {
    sum = _mm_add_pd(sum, _mm_mul_pd(
        _mm_loadu_pd(c + i),
        _mm_setr_pd(a[b[i]], a[b[i+1]])
    ));
}

if(n & 1) {
    sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1]));
}

double finalSum = _mm_cvtsd_f64(_mm_add_pd(
    sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1))
));

这会使用gcc -O2 -msse2 (4.4.1) 生成非常漂亮的汇编代码。

如您所知，具有偶数的n 将使此循环运行得更快以及对齐的c。如果您可以对齐c，请将_mm_loadu_pd 更改为_mm_load_pd 以获得更快的执行时间。

【讨论】：

不错，我忘了直接加载c了。
哦，嘿，哇——我不是故意要从@celion 那里抢走选定的答案...这是我个人的乐趣...
我相信我最终会克服它 :) 我认为我们两个答案的组合将是最佳的 - 我再次展开循环，并且您通过内在加载“c”。
通常您应该使用纯标量进行标量清理，例如与水平总和并行的finalSum = 0; 或finalsum = a[b[ n-1 ]] * c[ n-1 ];。

【解决方案2】：

我将从展开循环开始。类似的东西

double myNum1 = 0, myNum2=0;
for(int i=0;i<n;i+=2)
{
    myNum1 += a[b[ i ]] * c[ i ];
    myNum2 += a[b[i+1]] * c[i+1];
}
// ...extra code to handle the remainder when n isn't a multiple of 2...
double myNum = myNum1 + myNum2;

希望这允许编译器将负载与算术交错；配置文件并查看组件以查看是否有改进。理想情况下，编译器会生成 SSE 指令，但如果在实践中发生，我不会。

进一步展开可能会让您这样做：

__m128d sum0, sum1;
// ...initialize to zero...
for(int i=0;i<n;i+=4)
{
    double temp0 = a[b[ i ]] * c[ i ];
    double temp1 = a[b[i+1]] * c[i+1];
    double temp2 = a[b[i+2]] * c[i+2];
    double temp3 = a[b[i+3]] * c[i+3];
    __m128d pair0 = _mm_set_pd(temp0, temp1);
    __m128d pair1 = _mm_set_pd(temp2, temp3);
    sum0 = _mm_add_pd(sum0, pair0);
    sum1 = _mm_add_pd(sum1, pair1);
}
// ...extra code to handle the remainder when n isn't a multiple of 4...
// ...add sum0 and sum1, then add the result's components...

（对开头和结尾的伪代码表示歉意，我认为重要的部分是循环）。我不确定这是否会更快；这取决于各种延迟以及编译器重新排列所有内容的能力。确保您在之前和之后进行分析，看看是否有实际改进。

希望对您有所帮助。

【讨论】：

此外，它目前可能没有用，但我相信英特尔即将推出的架构 Larrabee 将具有收集/分散指令来处理这种情况。有关它的一些信息，请参阅drdobbs.com/architect/216402188?pgno=4。

【解决方案3】：

英特尔处理器可以发出两个浮点运算，但每个周期一个负载，因此内存访问是最严格的约束。考虑到这一点，我首先打算使用打包加载来减少加载指令的数量，并使用打包算术只是因为它方便。从那以后，我意识到内存带宽饱和可能是最大的问题，如果重点是让代码运行得更快而不是学习矢量化，那么所有与 SSE 指令有关的混乱可能都是过早的优化。

上海证券交易所

在不假设b 中的索引的情况下，尽可能少的负载需要展开循环四次。一个 128 位加载从 b 获得四个索引，两个 128 位加载每个从 c 获得一对相邻的双精度数，收集 a 需要独立的 64 位加载。对于串行代码，每四次迭代需要 7 个周期。（如果对a 的访问不能很好地缓存，足以使我的内存带宽饱和）。我省略了一些烦人的事情，比如处理不是 4 的倍数的迭代次数。

entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c)
  xorpd xmm0, xmm0
  xor r8, r8
loop:
  movdqa xmm1, [rdx+4*r8]
  movapd xmm2, [rcx+8*r8]
  movapd xmm3, [rcx+8*r8+8]
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm4, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm4, [rsi+8*r10]
  punpckhqdq xmm1, xmm1
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm5, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm5, [rsi+8*r10]
  add    r8,   4
  cmp    r8,   rdi
  mulpd  xmm2, xmm4
  mulpd  xmm3, xmm5
  addpd  xmm0, xmm2
  addpd  xmm0, xmm3
  jl loop

取出索引是最复杂的部分。 movdqa 从一个 16 字节对齐的地址加载 128 位整数数据（Nehalem 会因混合“整数”和“浮点”SSE 指令而产生延迟惩罚）。 punpckhqdq 将高 64 位移动到低 64 位，但在整数模式下，与更简单的名称 movhlpd 不同。 32 位移位在通用寄存器中完成。 movhpd 将一个 double 加载到 xmm 寄存器的上部而不影响下部 - 这用于将 a 的元素直接加载到打包寄存器中。

这段代码明显比上面的代码快，而上面的代码又比简单的代码快，而且在每个访问模式上，除了简单的情况B[i] = i，其中天真的循环实际上是最快的。我还尝试了一些东西，比如 Fortran 中 SUM(A(B(:)),C(:)) 周围的函数，它最终基本上等同于简单循环。

我在具有 4GB DDR2-667 内存、4 个模块的 Q6600（2.4Ghz 的 65 nm Core 2）上进行了测试。测试内存带宽大约为 5333 MB/s，所以看起来我只看到一个通道。我正在使用 Debian 的 gcc 4.3.2-1.1、-O3 -Ffast-math -msse2 -Ftree-vectorize -std=gnu99 进行编译。

为了测试，我让n 为一百万，初始化数组，使a[b[i]] 和c[i] 都等于1.0/(i+1)，并使用几种不同的索引模式。一个为a分配一百万个元素并将b设置为随机排列，另一个为a分配10M个元素并每10个使用一次，最后一个为a分配10M个元素并通过添加设置b[i+1]一个从 1 到 9 到 b[i] 的随机数。我正在计时使用gettimeofday 调用一次需要多长时间，通过在数组上调用clflush 来清除缓存，并测量每个函数的1000 次试验。我使用来自criterion 的一些代码（特别是statistics 包中的内核密度估计器）绘制了平滑的运行时分布。

带宽

现在，关于带宽的实际重要说明。 5333MB/s 和 2.4Ghz 时钟刚好超过每个周期两个字节。我的数据足够长，任何东西都不应该被缓存，如果一切都未命中，将循环的运行时间乘以每次迭代加载的 (16+2*16+4*64) 字节，几乎可以得到我的系统拥有的 ~5333MB/s 带宽.如果没有 SSE，应该很容易使带宽饱和。即使假设a 被完全缓存，只需读取b 和c 进行一次迭代就会移动12 个字节的数据，并且天真的可以通过流水线在第三个周期开始新的迭代。

假设在a 上没有完全缓存，那么算术和指令计数就不会成为瓶颈。如果我的代码中的大部分加速来自于向b 和c 发出更少的负载，我不会感到惊讶，因此可以有更多空间来跟踪和推测a 上过去的缓存未命中。

更广泛的硬件可能会产生更大的影响。运行三个 DDR3-1333 通道的 Nehalem 系统需要每个周期移动 3*10667/2.66 = 12.6 字节以使内存带宽饱和。如果a 适合缓存，这对于单个线程来说是不可能的 - 但是在 64 字节时，向量上的行缓存未命中会迅速加起来 - 缓存中缺少的循环中的四个负载中的一个将平均所需带宽提高到16 字节/周期。

【讨论】：

Sandybridge 及更高版本每个时钟可以执行 2 次加载。 Core 2 / Nehalem 也与后来的 Intel CPU 不同，movd/q reg, xmm 具有 3/clock（0.33c 吞吐量）。后来的 CPU（从 ~2013 年的 Haswell 开始）只有 1/clock 的吞吐量，所以那里会出现瓶颈。 agner.org/optimize/instruction_tables.pdf （所以你想用 mov / shr 进行 64 位标量加载）。但是，是的，这可能对 Core2 / Nehalem 有好处。
[rcx+8*r8+8] - 应该是+16，否则会错位并与之前的movapd重叠。

【解决方案4】：

简短的回答没有。长答案是的，但效率不高。您将因执行不对齐的负载而受到惩罚，这将抵消任何好处。除非您可以保证 b[i] 连续索引是对齐的，否则在矢量化之后您很可能会遇到更差的性能

如果您事先知道索引是什么，最好的办法是展开并指定显式索引。我使用模板专业化和代码生成做了类似的事情。如果你有兴趣，我可以分享

要回答您的评论，您基本上必须专注于数组。立即尝试的最简单的方法是将循环阻止两倍，分别加载低和高 a，然后像通常一样使用 mm*_pd。伪代码：

__m128d a, result;
for(i = 0; i < n; i +=2) {
  ((double*)(&a))[0] = A[B[i]];
  ((double*)(&a))[1] = A[B[i+1]];
  // you may also load B using packed integer instruction
  result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i])));
}

我不记得确切的函数名称，可能需要仔细检查。此外，如果您知道不存在别名问题，请在指针中使用 restrict 关键字。这将使编译器更具侵略性。

【讨论】：

你能解释一下我将如何进行矢量化（即使是不对齐的惩罚）？我很想自己对性能进行基准测试。
这行不通，因为索引的双重间接性。
我认为 restrict 在这里没有任何好处，因为所有写入都是针对局部变量的。如果他们计算 d[i] = a[b[i]] * c[i]，那么对 d 进行限制会有所帮助。

【解决方案5】：

由于数组索引的双重间接性，这不会按原样进行矢量化。由于您正在使用双打，因此从 SSE 中获得的收益很少或根本没有，特别是因为大多数现代 CPU 无论如何都有 2 个 FPU。

【讨论】：

错了，SSE2 允许在一个 128 位 SSE 寄存器中同时使用两个 64 位双精度。
@Liranuna - 在一个 SSE 寄存器中处理两个双精度如何比使用两个 FPU 更有优势？实际上，将两个不连续的双精度值打包到 SSE 寄存器中的额外开销几乎肯定会使 SSE 解决方案比标量实现慢。
@Paul：SSE 并不是一个神奇的优化护罩。如果你滥用它，你会发现它比幼稚的代码慢。然而，正确使用 SSE 总能让您的速度提升至少 30%。
@LiraNuna - 我知道 - 在这种情况下，我主张反对 SSE。
@Paul，为什么不进一步利用已经在使用的功能呢？ x86_64 的 fpu 已经与 mulsd 和 addsd 一致，为什么不使用相同的指令（花费相同数量的周期）来完成 DOUBLE 工作？请参阅我对如何在这种情况下充分利用 SSE(2) 的回答。