代码之家 › 专栏 › 技术社区 › Mike

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

vectorization simd sse x86-64 x86

Mike · 技术社区 · 16 年前

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
}

5 回复 | 直到 16 年前

phuclv 6 年前

这是我的尝试,经过充分优化和测试:

#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

n 将使此循环运行得更快,并且对齐 c . 如果你能对齐 C _mm_loadu_pd 到 _mm_load_pd 以获得更快的执行时间。

phuclv 6 年前

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...

希望有帮助。

Brandon 16 年前

上海证券交易所

尽可能少的荷载,且不假设 b 需要展开循环四次。一个128位加载从中获取四个索引 B ,两个128位加载,每个加载从 c ,及聚会 a 所需的独立64位加载。这是串行代码每四次迭代7个周期的下限(足够饱和我的内存带宽,如果访问 A.

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(:))

测试内存带宽大约为5333MB/s,所以看起来我只看到了一个通道。我使用Debian的GCC4.3.2-1.1,-O3-Ffast math-msse2-Ftree vectorize-std=gnu99进行编译。

为了测试,我让 n a[b[i]] 和 c[i] 两者相等 1.0/(i+1) ,有几种不同的指数模式。分配 A. 有一百万个元素和集合对于一个随机排列,另一个分配有10万个元素,每10个使用一次,最后一个分配 A. 10米的元素和设置 b[i+1] 通过将1到9之间的随机数添加到 b[i] gettimeofday clflush criterion statistics 包装)。

带宽

A. B 和 C 因为一次迭代会移动12字节的数据,而天真的人可以通过流水线第三个周期开始新的迭代。

假设服务器上没有完全缓存 A. B 和 C 因此,更多的空间可以自由跟踪和推测过去的缓存未命中 A.

更广泛的硬件可能会产生更大的差异。运行三个DDR3-1333通道的Nehalem系统需要每周期移动3*10667/2.66=12.6字节,以饱和内存带宽。如果 A. 适合缓存-但在64字节时,向量上的行缓存未命中会快速累积-在我的循环中,缓存中丢失的四个负载中,只有一个会将平均所需带宽提高到16字节/周期。

Anycorn 16 年前

简短回答否。冗长回答是,但效率不高。您将因执行非对齐负载而受到惩罚,这将否定任何类型的好处。除非您能保证b[i]连续索引对齐,否则向量化后的性能很可能会更差

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

嗯 *_像往常一样。伪代码:

__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关键字。这将使编译器更具攻击性。

Paul R 16 年前

这不会像现在这样矢量化,因为数组索引具有双重间接性。由于您使用的是Double,因此从SSE中几乎得不到什么好处,尤其是在大多数现代CPU都有2个FPU的情况下。