logo标签列表

在x86_64上是否可以对“myNum += a[b[i]] * c[i]”进行矢量化？

x86x86-64ssesimdvectorization
11

如果可能的话，在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中的索引分布是什么？ - MSN
只是好奇，下面的建议是否有助于加速您的代码？ - celion
5个回答
8
这是我的尝试，完全优化过并经过测试：

以下是我所做的完全优化和测试：

#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，可以进一步提高执行速度。

1
不错，我忘了直接加载c。 - celion
哦，嘿，哇——我不是想从@celion那里夺走被选中的答案...这只是我的个人乐趣... - LiraNuna
我相信最终我会克服这个问题的 :) 我认为结合我们两个答案的方法会是最优解 - 我重新展开循环，而你通过内置函数加载“c”。 - celion
通常情况下，您应该使用纯量进行标量清理，例如finalSum = 0;或finalsum = a[b[n-1]] * c[n-1];与水平求和并行。 - Peter Cordes
4
我会从展开循环开始。类似这样：
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将具有gather/scatter指令来处理这种情况。请参阅http://www.drdobbs.com/architect/216402188?pgno=4以获取有关它的一些信息。 - celion
2

英特尔处理器每个周期可以发出两个浮点运算，但只能进行一次加载，因此内存访问是最紧迫的限制。因此，我首先使用打包加载来减少加载指令的数量，并且仅仅因为方便而使用打包算术。后来我意识到饱和内存带宽可能是最大的问题，如果旨在使代码快速运行而不是学习矢量化，则所有与SSE指令相关的操作可能都是过早的优化。

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之外的每种访问模式下都是如此，其中天真循环实际上是最快的。我还尝试了一些类似Fortran中SUM(A(B(:)),C(:))的函数，最终基本等同于简单循环。
我在Q6600（65 nm Core 2 at 2.4Ghz）上进行了测试，内存为4GB的DDR2-667，分为4个模块。测试内存带宽约为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并使用每十个元素，最后一个分配了一千万个元素的a并通过将1到9之间的随机数添加到b[i]来设置b[i+1]。我用gettimeofday计时一次调用需要多长时间，通过对数组调用clflush来清除缓存，并测量每个函数的1000个试验。我使用从criterion（特别是statistics包中的核密度估计器）获取的一些代码绘制平滑的运行时分布图。
现在，有关带宽的实际重要说明。2.4GHz时钟下的5333MB/s仅为每个周期超过两个字节。我的数据足够长，以至于没有任何可以缓存的内容，并且如果所有内容都未命中，则将我的循环运行时间乘以每次迭代加载的(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在movd/q reg，xmm中的吞吐量为3/clock（0.33c），也与后来的英特尔CPU不同。后来的CPU（从Haswell开始，约为2013年）在这方面只有1/clock的吞吐量，因此会出现瓶颈。https://www.agner.org/optimize/instruction_tables.pdf（因此，您需要执行64位标量加载以使用`mov` / shr进行解包）。但是，对于Core2 / Nehalem来说，这可能是不错的选择。 - Peter Cordes
[rcx+8*r8+8] - 应该是 +16，否则它会错位并重叠前面的 movapd。 - Peter Cordes
0

简短的答案是否定的。长一点的答案是可以，但是不高效。你会承担非对齐加载的惩罚，这将抵消任何类型的好处。除非你能保证 b[i] 连续的索引是对齐的，否则向量化后性能很可能会更差。

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

回答你的评论，你基本上要集中在一个数组上。立即尝试的最简单的方法是将循环按两个因素分块，分别加载低位和高位，然后像通常一样使用 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 关键字。 这将使编译器更加积极。

1
你能解释一下我该如何进行向量化处理吗（即使有非对齐惩罚）？我很想自己测试一下性能。 - Mike
这不会起作用，因为索引的双重间接性。 - Paul R
1
我认为在这里使用restrict并没有任何好处，因为所有的写操作都是针对本地变量的。如果它们计算的是d[i] = a[b[i]] * c[i]，那么将d指定为restrict会有帮助。 - celion
0

由于数组索引的双重间接性，这样做不会进行向量化。由于你正在使用双精度浮点数，因此几乎没有从SSE中获得的好处，特别是现代大多数CPU都有2个FPU。

错误，SSE2允许在一个128位的SSE寄存器中同时处理两个64位的双精度浮点数。 - LiraNuna
@Liranuna - 如何在一个SSE寄存器中处理两个双精度浮点数比使用两个FPU更有优势？实际上，将两个不连续的双精度浮点数打包到一个SSE寄存器中的额外开销几乎肯定会使SSE解决方案比标量实现更慢。 - Paul R
@Paul：SSE并不是一个神奇的优化保护罩。如果你滥用它，你会发现它比朴素的代码更慢。然而，适当使用SSE将始终为您提供至少30%的速度提升。 - LiraNuna
@LiraNuna - 我知道 - 在这种情况下，我是反对使用SSE的。 - Paul R
@Liranuna - 你可能会看到一些小的好处，但我怀疑在这种情况下它不会接近2倍，因为算术操作的数量很少，并且每次迭代有很多负载。你还有一个未对齐的负载，除了Core i7之外的任何东西都很昂贵。不过，值得进行基准测试，并将吞吐量与简单的标量实现进行比较。 - Paul R
显示剩余3条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接