c++ - Add+Mul 使用 Intrinsics 变得更慢——我哪里错了？

Question

有这个数组:

alignas(16) double c[voiceSize][blockSize];

这是我要优化的功能:

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

    for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
        pC[sampleIndex] = value + deltaValue * sampleIndex;
    }
}

这是我的内在函数 (SSE2) 尝试:

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

    __m128d value_add = _mm_set1_pd(value);
    __m128d deltaValue_mul = _mm_set1_pd(deltaValue);

    for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex += 2) {
        __m128d result_mul = _mm_setr_pd(sampleIndex, sampleIndex + 1);
        result_mul = _mm_mul_pd(result_mul, deltaValue_mul);
        result_mul = _mm_add_pd(result_mul, value_add);

        _mm_store_pd(pC + sampleIndex, result_mul);
    }   
}

不幸的是，它比“标量”(即使是自动优化的)原始代码慢:)

您认为瓶颈在哪里？我哪里错了？

我正在使用 MSVC、Release/x86、/02 优化标志(支持快速代码) .

编辑:这样做(@wim 建议)，性能似乎比 C 版本更好:

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

    __m128d value_add = _mm_set1_pd(value);
    __m128d deltaValue_mul = _mm_set1_pd(deltaValue);

    __m128d sampleIndex_acc = _mm_set_pd(-1.0, -2.0);
    __m128d sampleIndex_add = _mm_set1_pd(2.0);

    for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex += 2) {
        sampleIndex_acc = _mm_add_pd(sampleIndex_acc, sampleIndex_add);
        __m128d result_mul = _mm_mul_pd(sampleIndex_acc, deltaValue_mul);
        result_mul = _mm_add_pd(result_mul, value_add);

        _mm_store_pd(pC + sampleIndex, result_mul);
    }
}

为什么？ _mm_setr_pd 贵吗？

Answer 1

Why? Is _mm_setr_pd expensive?

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

This is the function I'm trying to optimize:

您只是用线性函数填充数组。您每次在循环内都在重新乘以。这避免了对除整数循环计数器以外的任何事物的循环携带依赖性，但您会因在循环内执行大量工作而遇到吞吐量瓶颈。

即你在计算 a[i] = c + i*scale分别为每一步。但是，您可以将其强度降低为 a[i+n] = a[i] + (n*scale) .所以你只有一个 addpd每个结果 vector 的指令。

这将引入一些舍入误差，这些误差会累积而不是从头开始重做计算，但是 double无论如何，对于您正在做的事情来说可能有点矫枉过正。

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

因此，您需要展开多个 vector 以隐藏 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 .这可能与非常小的 block 大小相关，特别是如果您通过仅根据需要在读取该数组之前生成该数组的小块来缓存 block 代码。

无论如何，这会通过 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 addpd飞行中的指示。这对 Haswell 来说绰绰有余，但只有 Skylake 可以承受的一半。

仍然，存储吞吐量为每时钟 1 个 vector ，超过 1 addpd每个时钟没有用。 理论上，这可以以每个时钟周期大约 16 字节的速度运行，并使存储吞吐量饱和。即 1 vector/2 double每时钟 s。

具有更宽 vector (4 double s) 的 AVX 在 Haswell 及更高版本上仍可以每时钟 1 个 vector ，即每时钟 32 字节。 (假设输出数组在 L1d 缓存或什至 L2 缓存中是热的。)

---

更好的是:根本不要将这些数据存储在内存中；即时重新生成。

在需要时即时生成它，如果使用它的代码只读取它几次，并且也被手动矢量化。