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
}
4

5 回答 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)生成了非常漂亮的汇编代码。

正如你所知道的，有一个 evenn会使这个循环和对齐的一样快c。如果您可以 align c，请更改_mm_loadu_pd为_mm_load_pd以获得更快的执行时间。

于 2010-03-02T01:54:03.710 回答
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...

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

希望有帮助。

于 2010-02-28T11:28:16.880 回答
2

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

上证所

在不假设 in 索引的情况下，尽可能少的负载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，天真的循环实际上是最快的。我还尝试了一些东西，比如SUM(A(B(:)),C(:))Fortran 中的一个函数，它最终基本上等同于简单循环。

我在具有 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一百万，用一些不同的索引模式初始化数组 soa[b[i]]和c[i]both equal 。1.0/(i+1)一个分配a一百万个元素并设置b为随机排列，另一个分配a10M 个元素并每 10 个使用一次，最后一个分配a10M 个元素并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 字节/循环。

于 2010-03-05T13:21:46.837 回答
0

简短的回答没有。长答案是的，但效率不高。您将因执行不对齐的负载而受到惩罚，这将抵消任何好处。除非您可以保证 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 关键字。这将使编译器更具侵略性。

于 2010-02-28T05:00:40.740 回答
0

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

于 2010-02-28T09:14:16.050 回答