在 GCC 中生成没有 cmp 指令的循环

Question

我尝试使用 GCC 和内在函数优化许多紧密循环。例如下面的函数。

void triad(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    int i;
    __m256 k4 = _mm256_set1_ps(k);
    for(i=0; i<n; i+=8) {
        _mm256_store_ps(&z[i], _mm256_add_ps(_mm256_load_ps(&x[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i]))));
    }
}

这会产生一个像这样的主循环

20: vmulps ymm0,ymm1,[rsi+rax*1]
25: vaddps ymm0,ymm0,[rdi+rax*1]
2a: vmovaps [rdx+rax*1],ymm0
2f: add    rax,0x20
33: cmp    rax,rcx
36: jne    20 

但是cmp 指令是不必要的。我们可以将基指针（rsi、rdi 和 rdx）设置到数组的末尾并将 rax 设置为 -sizeof(float)*n，而不是让 rax 从零开始并在 sizeof(float)*n 结束然后测试为零。我可以用我自己的汇编代码来做到这一点

.L2  vmulps          ymm1, ymm2, [rdi+rax]
     vaddps          ymm0, ymm1, [rsi+rax]
     vmovaps         [rdx+rax], ymm0
     add             rax, 32
     jne             .L2

但我无法让 GCC 做到这一点。我现在有几个测试，这有很大的不同。直到最近 GCC 和内在函数已经很好地切断了我，所以我想知道是否有编译器开关或重新排序/更改我的代码的方法，因此 cmp 指令不是由 GCC 生成的。

我尝试了以下方法，但它仍然产生cmp。我尝试过的所有变体仍然产生cmp。

void triad2(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    float *x2 = x+n;
    float *y2 = y+n;
    float *z2 = z+n;    
    int i;
    __m256 k4 = _mm256_set1_ps(k);
    for(i=-n; i<0; i+=8) {
        _mm256_store_ps(&z2[i], _mm256_add_ps(_mm256_load_ps(&x2[i]), _mm256_mul_ps(k4, _mm256_load_ps(&y2[i]))));
    }
}

编辑： 我有兴趣为适合 L1 缓存的数组（实际上是 n=2048）最大化这些函数的指令级并行性 (ILP)。虽然展开可用于提高带宽，但它会降低 ILP（假设无需展开即可获得全部带宽）。

编辑： 这是 Core2（Nehalem 之前）、IvyBridge 和 Haswell 系统的结果表。 Intrinsics 是使用intrinsic 的结果，unroll1 是我没有使用cmp 的汇编代码，unroll16 是我的汇编代码展开16 次。百分比是峰值性能的百分比（frequency*num_bytes_cycle，其中 num_bytes_cycle 对于 SSE 是 24，对于 AVX 是 48，对于 FMA 是 96）。

                 SSE         AVX         FMA
intrinsic      71.3%       90.9%       53.6%      
unroll1        97.0%       96.1%       63.5%
unroll16       98.6%       90.4%       93.6%
ScottD         96.5%
32B code align             95.5%

对于 SSE，我在不展开的情况下得到的结果几乎与展开一样好，但前提是我不使用 cmp。在 AVX 上，我在不展开和不使用 cmp 的情况下获得最佳结果。有趣的是，在 IB 上展开实际上更糟。在 Haswell 上，我通过展开获得了迄今为止最好的结果。这就是为什么我问这个question。可以在那个问题中找到测试它的源代码。

编辑：

根据 ScottD 的回答，现在我的 Core2 系统（Nehalem 64 位之前的模式）的内部函数几乎达到了 97%。 我不确定为什么 cmp 确实很重要，因为它应该无论如何，每次迭代需要 2 个时钟周期。对于 Sandy Bridge，事实证明效率损失是由于代码对齐而不是额外的cmp。无论如何，在 Haswell 上只能展开工作。

Answer 1

Intel Ivy Bridge 或更高版本上的指令解码器可以将 cmp 和 jne 融合为流水线中的单个操作（称为宏操作融合），因此在这些最近的处理器上 cmp 无论如何都应该消失。

Answer 2

最终代码：

#define SF sizeof(float)
#ifndef NO                   //floats per vector, compile with -DNO = 1,2,4,8,...
#define NO 8                 //MUST be power of two
#endif

void triadfinaler(float const *restrict x, float const *restrict y,   \
                  float *restrict z, size_t n)
{
  float *restrict d = __builtin_assume_aligned(z, NO*SF);       //gcc builtin,
  float const *restrict m = __builtin_assume_aligned(y, NO*SF); //optional but produces
  float const *restrict a = __builtin_assume_aligned(x, NO*SF); //better code
  float const k = 3.14159f;
  n*=SF;
  while (n &= ~((size_t)(NO*SF)-1))    //this is why NO*SF must be power of two
    {
      size_t nl = n/SF;
      for (size_t i = 0; i<NO; i++)
        {
          d[nl-NO+i] = k * m[nl-NO+i] + a[nl-NO+i];
        }
      n -= (NO*SF);
    }
}

我更喜欢让编译器选择指令，而不是使用内在函数（尤其是因为您使用了 intel-intrinsics，而 gcc 并不喜欢它）。无论如何，以下代码在 gcc 4.8 上为我生成了不错的程序集：

void triad(float *restrict x, float *restrict y, float *restrict z, size_t n)
//I hope you weren't aliasing any function arguments... Oh, an it's void, not float
{
  float *restrict d = __builtin_assume_aligned(z, 32);  // Uh, make sure your arrays
  float *restrict m = __builtin_assume_aligned(y, 32);  // are aligned? Faster that way
  float *restrict a = __builtin_assume_aligned(x, 32);  //
  float const k = 3.14159f;
  while (n &= ~((size_t)0x7))       //black magic, causes gcc to omit code for non-multiples of 8 floats
    {
      n -= 8;                       //You were always computing on 8 floats at a time, right?
      d[n+0] = k * m[n+0] + a[n+0]; //manual unrolling
      d[n+1] = k * m[n+1] + a[n+1];
      d[n+2] = k * m[n+2] + a[n+2];
      d[n+3] = k * m[n+3] + a[n+3];
      d[n+4] = k * m[n+4] + a[n+4];
      d[n+5] = k * m[n+5] + a[n+5];
      d[n+6] = k * m[n+6] + a[n+6];
      d[n+7] = k * m[n+7] + a[n+7];
    }
}

这为我的 corei7avx2 生成了不错的代码，带有 -O3：

triad:
    andq    $-8, %rcx
    je  .L8
    vmovaps .LC0(%rip), %ymm1

.L4:
    subq    $8, %rcx
    vmovaps (%rsi,%rcx,4), %ymm0
    vfmadd213ps (%rdi,%rcx,4), %ymm1, %ymm0
    vmovaps %ymm0, (%rdx,%rcx,4)
    andq    $-8, %rcx
    jne .L4
    vzeroupper
.L8:
    rep ret
    .cfi_endproc

.LC0:
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000

编辑： 我对编译器没有将这段代码优化到最后一条指令感到有点失望，所以我把它弄得更乱了。只需更改循环中的事物顺序就可以摆脱编译器发出的AND，这让我走上了正确的轨道。然后我只需要让它不在循环中进行不必要的地址计算。叹息。

void triadtwo(float *restrict x, float *restrict y, float *restrict z, size_t n)
{
  float *restrict d = __builtin_assume_aligned(z, 32);
  float *restrict m = __builtin_assume_aligned(y, 32);
  float *restrict a = __builtin_assume_aligned(x, 32);
  float const k = 3.14159f;
  n<<=2;
  while (n &= -32)
    {
      d[(n>>2)-8] = k * m[(n>>2)-8] + a[(n>>2)-8];
      d[(n>>2)-7] = k * m[(n>>2)-7] + a[(n>>2)-7];
      d[(n>>2)-6] = k * m[(n>>2)-6] + a[(n>>2)-6];
      d[(n>>2)-5] = k * m[(n>>2)-5] + a[(n>>2)-5];
      d[(n>>2)-4] = k * m[(n>>2)-4] + a[(n>>2)-4];
      d[(n>>2)-3] = k * m[(n>>2)-3] + a[(n>>2)-3];
      d[(n>>2)-2] = k * m[(n>>2)-2] + a[(n>>2)-2];
      d[(n>>2)-1] = k * m[(n>>2)-1] + a[(n>>2)-1];
      n -= 32;
    }
}

丑陋的代码？是的。但是大会：

triadtwo:
    salq    $2, %rcx
    andq    $-32, %rcx
    je  .L54
    vmovaps .LC0(%rip), %ymm1

.L50:
    vmovaps -32(%rsi,%rcx), %ymm0
    vfmadd213ps -32(%rdi,%rcx), %ymm1, %ymm0
    vmovaps %ymm0, -32(%rdx,%rcx)
    subq    $32, %rcx
    jne .L50
    vzeroupper
.L54:
    rep ret
    .cfi_endproc
.LC0:
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000
    .long   1078530000

Mmmmhhh，光荣的循环中的五个指令，宏操作可熔减和分支......

Answer 3

这个怎么样。编译器是 gcc 4.9.0 mingw x64：

void triad(float *x, float *y, float *z, const int n) {
    float k = 3.14159f;
    intptr_t i;
    __m256 k4 = _mm256_set1_ps(k);

    for(i = -n; i < 0; i += 8) {
        _mm256_store_ps(&z[i+n], _mm256_add_ps(_mm256_load_ps(&x[i+n]), _mm256_mul_ps(k4, _mm256_load_ps(&y[i+n]))));
    }
}

gcc -c -O3 -march=corei7 -mavx2 triad.c

0000000000000000 <triad>:
   0:   44 89 c8                mov    eax,r9d
   3:   f7 d8                   neg    eax
   5:   48 98                   cdqe
   7:   48 85 c0                test   rax,rax
   a:   79 31                   jns    3d <triad+0x3d>
   c:   c5 fc 28 0d 00 00 00 00 vmovaps ymm1,YMMWORD PTR [rip+0x0]
  14:   4d 63 c9                movsxd r9,r9d
  17:   49 c1 e1 02             shl    r9,0x2
  1b:   4c 01 ca                add    rdx,r9
  1e:   4c 01 c9                add    rcx,r9
  21:   4d 01 c8                add    r8,r9

  24:   c5 f4 59 04 82          vmulps ymm0,ymm1,YMMWORD PTR [rdx+rax*4]
  29:   c5 fc 58 04 81          vaddps ymm0,ymm0,YMMWORD PTR [rcx+rax*4]
  2e:   c4 c1 7c 29 04 80       vmovaps YMMWORD PTR [r8+rax*4],ymm0
  34:   48 83 c0 08             add    rax,0x8
  38:   78 ea                   js     24 <triad+0x24>

  3a:   c5 f8 77                vzeroupper
  3d:   c3                      ret

就像您的手写代码一样，gcc 使用 5 条指令进行循环。 gcc 代码使用 scale=4，而您的代码使用 scale=1。我能够让 gcc 在 5 指令循环中使用 scale=1，但是 C 代码很笨拙，循环中的 2 条 AVX 指令从 5 字节增长到 6 字节。