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

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

（在开始和结束时为伪代码道歉，我认为重要的部分是循环）。我不确定这是否会更快;它取决于各种延迟以及编译器重新排列所有内容的程度。确保您在前后查看是否有实际改进。

希望有所帮助。

英特尔处理器可以发出两个浮点运算，但每个周期一次加载，因此内存访问是最严格的约束。考虑到这一点，我首先瞄准使用打包负载来减少加载指令的数量，并且使用打包算法只是因为它很方便。我已经意识到饱和内存带宽可能是最大的问题，如果关键是要使代码快速而不是学习矢量化，那么所有乱用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天真的循环实际上是最快的。我也尝试了一些类似于函数的东西 SUM(A(B(:)),C(:)) 在Fortran中，它最终基本上等同于简单的循环。

我测试了Q6600（65 nm Core 2，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 10M元素，每10次使用，最后一次分配 a 拥有10M元素和设置 b[i+1] 通过添加从1到9的随机数 b[i]。我计时一次通话需要多长时间 gettimeofday，通过调用清除缓存 clflush 在阵列上，并测量每个功能的1000个试验。我使用了一些代码来绘制平滑的运行时分布 标准 （特别是，核心密度估计 statistics 包）。

带宽

现在，关于带宽的实际重要说明。具有2.4Ghz时钟的5333MB / s每个周期仅超过两个字节。我的数据足够长，没有任何东西可以缓存，并且如果一切都没有，我的循环的运行时间乘以每次迭代加载的（16 + 2 * 16 + 4 * 64）个字节，这几乎可以让我的系统具有~5333MB / s的带宽。 。没有SSE就可以很容易地使带宽饱和。即使假设 a 完全缓存，只是阅读 b 和 c 对于一次迭代，移动12个字节的数据，并且naive可以使用流水线技术从第三个循环开始新的迭代。

假设任何不完整的缓存 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])));
}

我不记得确切的功能名称，可能要仔细检查。 另外，如果您知道不存在别名问题，请使用带有指针的restrict关键字。 这将使编译器更具攻击性。

由于数组索引的双重间接，因此不会按原样进行矢量化。由于你正在使用双打，因此SSE几乎没有任何东西可以获得，特别是因为大多数现代CPU都有2个FPU。