add+mul在使用intrinsics时变慢了-我错在哪里?

为什么?是不是很贵?

有点;至少需要一次洗牌。更重要的是,在这种情况下,计算每个标量操作数是昂贵的,而且正如@spectras'的答案所示,gcc至少无法将其自动矢量化为 paddd / cvtdq2pd . 相反,它从标量整数重新计算每个操作数,执行 int -gt; double 分别转换,然后将它们混合在一起。

这是我试图优化的功能:

你只是用一个线性函数填充一个数组 . 你在循环中每次都在增加。这避免了循环对整数循环计数器以外的任何东西的依赖性,但是在循环中做这么多工作会遇到吞吐量瓶颈。

也就是说,你在计算 a[i] = c + i*scale 每一步都要分开。而是相反 你可以把力量降低到 a[i+n] = a[i] + (n*scale) . 所以你只有一个 addpd 每个结果向量的指令。

这将引入一些舍入误差,这些误差累积起来,而不是从头开始重新计算,但是 双重的 可能对你的所作所为有点过分了。

它还以引入对fp add而不是整数的串行依赖为代价。但是您已经在“优化”版本中有一个带循环的fp-add依赖链,它使用 sampleIndex_acc = _mm_add_pd(sampleIndex_acc, sampleIndex_add); 在循环内,使用fp+=2.0而不是从整数重新转换。

因此,您需要使用多个向量展开以隐藏fp延迟,并同时保持至少3或4个fp添加。 . (Haswell:3周期延迟,每时钟一次吞吐量。Skylake:4周期延迟,每时钟吞吐量2。)另请参见 Why does mulss take only 3 cycles on Haswell, different from Agner's instruction tables? 有关使用多个累加器展开循环执行依赖项(点积)的类似问题的详细信息。

void Process(int voiceIndex, int blockSize) {    
    double *pC = c[voiceIndex];
    double val0 = start + step * delta;
    double deltaValue = rate * delta;

    __m128d vdelta2 = _mm_set1_pd(2 * deltaValue);
    __m128d vdelta4 = _mm_add_pd(vdelta2, vdelta2);

    __m128d v0 = _mm_setr_pd(val0, val0 + deltaValue);
    __m128d v1 = _mm_add_pd(v0, vdelta2);
    __m128d v2 = _mm_add_pd(v0, vdelta4);
    __m128d v3 = _mm_add_pd(v1, vdelta4);

    __m128d vdelta8 = _mm_mul_pd(vdelta2, _mm_set1_pd(4.0));

    double *endp = pC + blocksize - 7;  // stop if there's only room for 7 or fewer doubles
      // or use -8 and have your cleanup handle lengths of 1..8
      // since the inner loop always calculates results for next iteration
    for (; pC < endp ; pC += 8) {
        _mm_store_pd(pC, v0);
        v0 = _mm_add_pd(v0, vdelta8);

        _mm_store_pd(pC+2, v1);
        v1 = _mm_add_pd(v1, vdelta8);

        _mm_store_pd(pC+4, v2);
        v2 = _mm_add_pd(v2, vdelta8);

        _mm_store_pd(pC+6, v3);
        v3 = _mm_add_pd(v3, vdelta8);
    }
    // if (blocksize % 8 != 0) ... store final vectors
}

建立时是加还是乘的选择 vdelta4 / vdelta8 并不是非常重要;我试图避免太久的依赖链之前的第一个商店可能发生。自从 v0 通过 v3 也需要计算一下,创建一个 VDelTA4 而不只是做一个链子 v2 = v1+vdelta2 . 也许最好是创造 VDelTA4 乘从 4.0*delta ,然后加倍 VDelTA8 . 这可能与非常小的块大小有关,特别是如果在读取代码之前仅根据需要生成此数组的小块来缓存块代码。

无论如何,这将编译成一个非常有效的gcc和msvc内部循环( on the Godbolt compiler explorer )

;; MSVC -O2
$LL4@Process:                    ; do {
    movups  XMMWORD PTR [rax], xmm5
    movups  XMMWORD PTR [rax+16], xmm0
    movups  XMMWORD PTR [rax+32], xmm1
    movups  XMMWORD PTR [rax+48], xmm2
    add     rax, 64                             ; 00000040H
    addpd   xmm5, xmm3              ; v0 += vdelta8
    addpd   xmm0, xmm3              ; v1 += vdelta8
    addpd   xmm1, xmm3              ; v2 += vdelta8
    addpd   xmm2, xmm3              ; v3 += vdelta8
    cmp     rax, rcx
    jb      SHORT $LL4@Process   ; }while(pC < endp)

它有4个独立的依赖链,通过xmm0、1、2和5。所以有足够的指令级并行度来保持4 ADPD 飞行指令。这对Haswell来说已经绰绰有余了,但Skylake可以维持的一半。

不过,如果存储吞吐量为每个时钟1个矢量,则大于1 ADPD 每个钟都没用。 理论上,每时钟周期可以运行大约16字节,并且饱和存储吞吐量。 即1矢量/2 双重的 时钟。

矢量更宽的AVX(4 双重的 s)仍然可以在haswell和更高版本上以每个时钟1个矢量运行,即每个时钟32字节。(假设输出数组在L1D缓存中是热的,甚至可能是二级缓存。)

更好的方法是:根本不要将这些数据存储在内存中;动态地重新生成。

在需要时动态生成它,如果使用它的代码只读取它几次,而且也是手动矢量化的。

在我的系统上, g++ test.cpp -march=native -O2 -c -o test

这将输出正常版本(循环体提取):

  30:   c5 f9 57 c0             vxorpd %xmm0,%xmm0,%xmm0
  34:   c5 fb 2a c0             vcvtsi2sd %eax,%xmm0,%xmm0
  38:   c4 e2 f1 99 c2          vfmadd132sd %xmm2,%xmm1,%xmm0
  3d:   c5 fb 11 04 c2          vmovsd %xmm0,(%rdx,%rax,8)
  42:   48 83 c0 01             add    $0x1,%rax
  46:   48 39 c8                cmp    %rcx,%rax
  49:   75 e5                   jne    30 <_Z11ProcessAutoii+0x30>

对于intrinsics版本:

  88:   c5 f9 57 c0             vxorpd %xmm0,%xmm0,%xmm0
  8c:   8d 50 01                lea    0x1(%rax),%edx
  8f:   c5 f1 57 c9             vxorpd %xmm1,%xmm1,%xmm1
  93:   c5 fb 2a c0             vcvtsi2sd %eax,%xmm0,%xmm0
  97:   c5 f3 2a ca             vcvtsi2sd %edx,%xmm1,%xmm1
  9b:   c5 f9 14 c1             vunpcklpd %xmm1,%xmm0,%xmm0
  9f:   c4 e2 e9 98 c3          vfmadd132pd %xmm3,%xmm2,%xmm0
  a4:   c5 f8 29 04 c1          vmovaps %xmm0,(%rcx,%rax,8)
  a9:   48 83 c0 02             add    $0x2,%rax
  ad:   48 39 f0                cmp    %rsi,%rax
  b0:   75 d6                   jne    88 <_Z11ProcessSSE2ii+0x38>

简而言之:编译器从C版本自动生成AVX代码。

在使用标志播放更多内容后进行编辑,以便仅在两种情况下使用SSE2:

g++ test.cpp -msse2 -O2 -c -o test

编译器仍然执行与使用内部函数生成的不同的操作。编译器版本:

  30:   66 0f ef c0             pxor   %xmm0,%xmm0
  34:   f2 0f 2a c0             cvtsi2sd %eax,%xmm0
  38:   f2 0f 59 c2             mulsd  %xmm2,%xmm0
  3c:   f2 0f 58 c1             addsd  %xmm1,%xmm0
  40:   f2 0f 11 04 c2          movsd  %xmm0,(%rdx,%rax,8)
  45:   48 83 c0 01             add    $0x1,%rax
  49:   48 39 c8                cmp    %rcx,%rax
  4c:   75 e2                   jne    30 <_Z11ProcessAutoii+0x30>

内部版本:

  88:   66 0f ef c0             pxor   %xmm0,%xmm0
  8c:   8d 50 01                lea    0x1(%rax),%edx
  8f:   66 0f ef c9             pxor   %xmm1,%xmm1
  93:   f2 0f 2a c0             cvtsi2sd %eax,%xmm0
  97:   f2 0f 2a ca             cvtsi2sd %edx,%xmm1
  9b:   66 0f 14 c1             unpcklpd %xmm1,%xmm0
  9f:   66 0f 59 c3             mulpd  %xmm3,%xmm0
  a3:   66 0f 58 c2             addpd  %xmm2,%xmm0
  a7:   0f 29 04 c1             movaps %xmm0,(%rcx,%rax,8)
  ab:   48 83 c0 02             add    $0x2,%rax
  af:   48 39 f0                cmp    %rsi,%rax
  b2:   75 d4                   jne    88 <_Z11ProcessSSE2ii+0x38>

编译器不在此展开循环。这可能是更好或更糟取决于许多事情。你可能想把这两个版本都试一下。