是否可以向量化 myNum += a[b[i]] * c[i];在 x86_64 上？

StackOverflow https://stackoverflow.com/questions/2350116

题

我将使用哪些内在函数在 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
}
有帮助吗？

解决方案

下面是我走在它，完全优化和测试：

#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为更快的执行时间。

其他提示

我将从展开循环开始。就像是

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 指令，但如果实际发生这种情况，我不会生成 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...

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

希望有帮助。

英特尔处理器可以发出两个浮点运算，但每个周期一次加载，因此内存访问是最严格的约束。考虑到这一点，我的目标首先是使用打包加载来减少加载指令的数量，并使用打包算术只是因为它很方便。从那以后我意识到内存带宽饱和可能是最大的问题，如果重点是让代码运行得更快而不是学习矢量化，那么所有对 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 将一个双精度数加载到 xmm 寄存器的上部，而不影响下部 - 这用于加载以下元素 a 直接进入打包寄存器。

此代码明显比上面的代码快，而上面的代码又比简单代码快，并且在除简单情况之外的每种访问模式上 B[i] = i 其中朴素循环实际上是最快的。我还尝试了一些类似功能的东西 SUM(A(B(:)),C(:)) 在 Fortran 中，最终基本上等同于简单循环。

我在具有 4GB DDR2-667 内存、4 个模块的 Q6600（65 nm Core 2，2.4Ghz）上进行了测试。测试内存带宽约为 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 次试验。我使用一些来自内部的代码绘制了平滑的运行时分布 标准 （特别是，核密度估计器 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 字节/周期。

短答案没有。龙回答是，但不是有效的。你会招致这样做，这将否定任何形式的利益不结盟负荷的惩罚。除非你能保证B [I]连续指数排列，你将极有可能在矢量

后表现更差

如果你预先知道指数，你最好是解开，并指定明确的指标。我没有类似的东西使用模板特殊化和代码生成。如果你有兴趣，我可以分享

回答您的评论，你基本上要集中一个阵列上。试试就最简单的办法是两个，负荷低，分别高一的一个因素阻止你循环，然后使用的毫米的* _ 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])));
}

我不记得函数的名称完全相同，可能要仔细检查。 此外，使用限制关键字与指针，如果你知道不可能有锯齿的问题。 这将允许编译器要更加积极。

这是不会进行向量化，因为该数组索引的双间接的，因为它是，。既然你用双打的工作很少有或没有从SSE获得，特别是最现代的CPU有2周的FPU反正。

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top