使用Intrinsics后,Add+Mul变慢了-我错在哪里?

4

以下是一个数组示例:

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);
    }   
}

很遗憾,即使自动优化,它比“scalar”代码慢。

在您的意见中,瓶颈在哪里?我哪里做错了?

我正在使用MSVCRelease/x86/02 优化标志(Favor fast code)。

编辑:按照@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 很耗费资源?

3
你的硬件支持AVX吗?那么优化过的代码可能会使用AVX指令集。 - Max
4
你能发布一些我们可以编译的内容吗?编译器可能会为你的平台生成更好的内部函数! - Matthieu Brucher
2
我会从类似于__m128d sampleIndex_vec = _mm_set_pd(-1.0,-2.0);__m128d sampleIndex_add = _mm_set1_pd(2.0);的东西开始,放在循环外面。在循环内部,您可以将__m128d result_mul = _mm_setr_pd(sampleIndex, sampleIndex + 1);替换为sampleIndex_vec = _mm_add_pd(sampleIndex_vec, sampleIndex_add);result_mul = sampleIndex_vec。这样就可以摆脱讨厌的_mm_setr_pd(sampleIndex, sampleIndex + 1);了。(未经测试。) - wim
2
@markzzz 我的第一条评论有误,我已将其删除并写了一条新的评论。(检查结果以查看您的SIMD代码是否正确)思路是您完全在SIMD域中更新sampleIndex计数器,这比每次迭代进行两个整数转换为双精度浮点数要高效得多。使用gcc -S查看两个版本之间的差异。 - wim
1
@markzzz 我不使用MSVC,但也许这个问题对你有用。我强烈建议在编写SIMD内嵌代码时检查生成的汇编代码。 - wim
显示剩余9条评论
2个回答

3

在我的系统上,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>

对于内在版本:

  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>

编译器不会在这里展开循环。这可能会有好坏之分,取决于许多因素。你可能想对两个版本进行基准测试。

1
是的,它会。使用“-march”选项指定目标架构。您甚至可以编译多个版本,并根据检测到的CPU功能在运行时选择其中一个。 - spectras
3
它可能表现出与预期不同的行为,但不一定会“崩溃”。编译器标志“-march=native”表示“将其编译为在与编译器相同的处理器上运行”。如果您有特定目标,请考虑精确编译到该目标。 - Caleth
2
你的汇编代码中没有使用AVX512,而是使用AVX1 VEX编码,带有2或3字节的前缀(第一个字节为“c4”或“c5”),而不是4字节的EVEX前缀。但是,自动向量化比手动向量化的代码更好,因为手动向量化会导致编译器在循环内部进行标量的int->fp转换,从而走入错误的路径。问题在于手动向量化非常糟糕。 - Peter Cordes
1
@markzzz:是的,众所周知,gcc优化比大多数情况下的MSVC更好。MSVC通常是4个主要x86编译器(clang,gcc,ICC,MSVC)中最糟糕的一款,而且它在制作针对某台机器优化的二进制文件时可选项非常有限。 - Peter Cordes
1
@markzzz:当然,可以使用2、3或4个sampleIndex_acc变量来隐藏FP-add延迟和瓶颈,从而提高吞吐量。即使用多个累加器展开。 - Peter Cordes
显示剩余10条评论

3

为什么 _mm_setr_pd 这么耗费资源?

有点耗费资源,至少需要进行一次洗牌操作。更重要的是,在这种情况下,计算每个标量操作都很耗费资源,并且正如 @spectras 的回答所示,至少在gcc中无法自动将其向量化为 paddd/cvtdq2pd。相反,它会从标量整数重新计算每个操作数,分别进行 int->double 转换,然后将它们混合在一起。

我正在尝试优化的函数:

你只是用线性函数填充一个数组。你每次都在循环内重新乘以倍数。这避免了除整数循环计数器之外的任何循环依赖,但由于在循环内做了这么多工作,你会遇到吞吐量瓶颈。

即你每次都在单独计算 a[i] = c + i*scale。但是,你可以通过强度减少将其简化为 a[i+n] = a[i] + (n*scale)。因此,每个结果向量只有一个 addpd 指令。

这将引入一些舍入误差,与重新从头计算相比会累积误差,但是对于你正在做的事情来说,double 可能过于复杂了。

它还带来了一个串行依赖 FP 加法的成本,而不是整数。但是你已经在你的“优化”版本中有一个循环 FP 加法依赖链,使用 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 时,选择加法或乘法并不是非常重要的;我试图避免太长的依赖链,在第一次存储之前就可以发生。由于也需要计算出 v0v3,因此创建一个 vdelta4 而不仅仅是创建一个 v2=v1+vdelta2 的链条似乎更有意义。也许更好的方法是用从 4.0*delta 开始的乘法来创建 vdelta4,然后将其加倍以获得 vdelta8。这对于非常小的块大小可能很重要,特别是如果您通过仅在需要之前生成此数组的小块并进行缓存块处理来优化您的代码,就在它将被读取之前。

无论如何,这在gcc和MSVC中都编译为非常高效的内部循环(在Godbolt编译器浏览器上)。

;; 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个向量,超过1个addpd每个时钟周期是没有用的。理论上,这可以以大约16字节每个时钟周期的速度运行,并饱和存储吞吐量。即每个时钟周期1个向量/2个doubles。
使用更宽的向量(4个doubles)的AVX仍然可以在Haswell及更高版本上以每个时钟周期1个向量的速度运行,即32字节每个时钟周期。(假设输出数组在L1d缓存或甚至L2中是热的。)

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

当需要时动态生成它,如果使用它的代码只读取几次,并且也手动进行了向量化。


2
很高兴看到有更有经验的人接替我所停留的地方,并且走得更远。这就是让我们一步步变得更好的事情。 - spectras
1
@markzzz:在这种情况下,您至少要为vdeltaX本身留下一个寄存器。如果您将所有8个用于累加器,则编译器必须溢出某些内容,或者是vdeltaX并执行addps xmm0,[esp] / addps xmm1,[esp] / ...(每个addps中的额外负载),或者是溢出/重新加载其中一个v0..7累加器,通过存储/重新加载延长循环传递依赖链,从而破坏了练习的目的。它不必是2的幂,但是:您可以使用6或7个累加器v0..v5和pC += 12展开。值得一试,可能会有所帮助。但是,然后您可能需要编写清理代码。 - Peter Cordes
1
当您定义alignas(16) double c[13][blocksize]时,如果您想使用125的块大小,则仍然像c[13][126]一样定义它,这样每行末尾都有一个未使用的元素。您可以定义一个宏,如#define roundup2(x) ((x + 1UL) & -2UL),这样您就可以执行alignas(16) double c[13][roundup2(blocksize)] - Peter Cordes
1
@markzzz:是的,我认为-7是对的。-8会留下一个完整的4向量集合而不会超出边缘。但根据您的清理代码,您实际上可能需要-8。内部循环计算下一次迭代的完整结果集,并且在最后一次迭代中您不需要它。 - Peter Cordes
1
循环退出分支通常会在最后一次迭代上预测错误(在那里它会掉落而不是被采取),因此有足够的时间运行那些没有等待结果的额外ADDPD指令。或者您可以向下舍入并执行一定数量的清理存储块,但是在现代CPU上,乱序执行会重叠执行所有操作,所以您当前的解决方案基本上很好。 - Peter Cordes
显示剩余13条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接