展开循环并使用矢量化进行独立求和答案

【问题标题】：Unroll loop and do independent sum with vectorization展开循环并使用矢量化进行独立求和
【发布时间】：2016-01-07 10:02:52
【问题描述】：

对于以下循环，如果我告诉它使用关联数学，GCC 只会对循环进行矢量化，例如-Ofast。

float sumf(float *x)
{
  x = (float*)__builtin_assume_aligned(x, 64);
  float sum = 0;
  for(int i=0; i<2048; i++) sum += x[i];
  return sum;
}

这是带有-Ofast -mavx的程序集

sumf(float*):
    vxorps  %xmm0, %xmm0, %xmm0
    leaq    8192(%rdi), %rax
.L2:
    vaddps  (%rdi), %ymm0, %ymm0
    addq    $32, %rdi
    cmpq    %rdi, %rax
    jne .L2
    vhaddps %ymm0, %ymm0, %ymm0
    vhaddps %ymm0, %ymm0, %ymm1
    vperm2f128  $1, %ymm1, %ymm1, %ymm0
    vaddps  %ymm1, %ymm0, %ymm0
    vzeroupper
    ret

这清楚地表明循环已被矢量化。

但是这个循环也有一个依赖链。为了克服加法的延迟，我需要在 x86_64 上展开并至少进行 3 次部分求和（不包括需要展开 8 次的 Skylake 以及需要在 Haswell 和 Broadwell 上展开 10 次的 FMA 指令进行加法） .据我了解，我可以使用-funroll-loops 展开循环。

这是带有-Ofast -mavx -funroll-loops 的程序集。

sumf(float*):
    vxorps  %xmm7, %xmm7, %xmm7
    leaq    8192(%rdi), %rax
.L2:
    vaddps  (%rdi), %ymm7, %ymm0
    addq    $256, %rdi
    vaddps  -224(%rdi), %ymm0, %ymm1
    vaddps  -192(%rdi), %ymm1, %ymm2
    vaddps  -160(%rdi), %ymm2, %ymm3
    vaddps  -128(%rdi), %ymm3, %ymm4
    vaddps  -96(%rdi), %ymm4, %ymm5
    vaddps  -64(%rdi), %ymm5, %ymm6
    vaddps  -32(%rdi), %ymm6, %ymm7
    cmpq    %rdi, %rax
    jne .L2
    vhaddps %ymm7, %ymm7, %ymm8
    vhaddps %ymm8, %ymm8, %ymm9
    vperm2f128  $1, %ymm9, %ymm9, %ymm10
    vaddps  %ymm9, %ymm10, %ymm0
    vzeroupper
    ret

GCC 确实展开了八次循环。但是，它不进行独立求和。它做了八个相关的和。这毫无意义，也比不展开更好。

如何让 GCC 展开循环并进行独立的部分求和？

编辑：

即使没有 SSE 的 -funroll-loops，Clang 也会展开为四个独立的部分总和，但我不确定它的 AVX 代码是否同样有效。无论如何，编译器不应该需要 -funroll-loops 和 -Ofast，所以很高兴看到 Clang 至少在 SSE 中这样做是正确的。

Clang 3.5.1 与 -Ofast。

sumf(float*):                              # @sumf(float*)
    xorps   %xmm0, %xmm0
    xorl    %eax, %eax
    xorps   %xmm1, %xmm1
.LBB0_1:                                # %vector.body
    movups  (%rdi,%rax,4), %xmm2
    movups  16(%rdi,%rax,4), %xmm3
    addps   %xmm0, %xmm2
    addps   %xmm1, %xmm3
    movups  32(%rdi,%rax,4), %xmm0
    movups  48(%rdi,%rax,4), %xmm1
    addps   %xmm2, %xmm0
    addps   %xmm3, %xmm1
    addq    $16, %rax
    cmpq    $2048, %rax             # imm = 0x800
    jne .LBB0_1
    addps   %xmm0, %xmm1
    movaps  %xmm1, %xmm2
    movhlps %xmm2, %xmm2            # xmm2 = xmm2[1,1]
    addps   %xmm1, %xmm2
    pshufd  $1, %xmm2, %xmm0        # xmm0 = xmm2[1,0,0,0]
    addps   %xmm2, %xmm0
    retq

带有-O3 的ICC 13.0.1 展开为两个独立的部分和。 ICC 显然只假设关联数学与 -O3。

.B1.8: 
    vaddps    (%rdi,%rdx,4), %ymm1, %ymm1                   #5.29
    vaddps    32(%rdi,%rdx,4), %ymm0, %ymm0                 #5.29
    vaddps    64(%rdi,%rdx,4), %ymm1, %ymm1                 #5.29
    vaddps    96(%rdi,%rdx,4), %ymm0, %ymm0                 #5.29
    addq      $32, %rdx                                     #5.3
    cmpq      %rax, %rdx                                    #5.3
    jb        ..B1.8        # Prob 99%                      #5.3

【问题讨论】：

手动添加8个累加器？
@user3528438，这违背了让编译器为我执行此操作的全部目的。无论如何我只会展开四次，如果我必须手动展开，我还不如使用内在函数（无论如何我都会在实践中这样做）。 ICC 偶然展开为两个部分总和。 ICC 更好。
我试过#pragma omp simd reduction(+:sum) aligned(x:64)和-fopenmp。那肯定做了更多的事情，但我无法阅读足够多的程序集来判断它是否解决了您的问题。可以吗？
公平地说，我认为您对编译器的要求太高了。也许再给它几年？
@Zboson：Skylake 有 4 个周期的延迟 vaddps，每个周期的吞吐量为 2。（它丢弃了 3c FP 加法单元，并使用 4 周期延迟 FMA 单元进行加法和乘法运算。）您需要 8 个矢量累加器来使 Skylake 的加法、mul 或 fma 的 FP 吞吐量饱和。我完全同意，如果编译器展开更聪明地使用更多累加器，那将是非常好的。 clang 3.7 on godbolt 使用 4，但毫无意义地展开更多。（uop 缓存很小，因此仅根据需要展开。gcc 默认仅使用-fprofile-use 展开。）

标签： c gcc x86 loop-unrolling auto-vectorization

【解决方案1】：

gcc 内部函数和__builtin_ 的一些使用会产生这样的结果：

typedef float v8sf __attribute__((vector_size(32)));
typedef uint32_t v8u32 __attribute__((vector_size(32)));

static v8sf sumfvhelper1(v8sf arr[4])
{
  v8sf retval = {0};
  for (size_t i = 0; i < 4; i++)
    retval += arr[i];
  return retval;
}

static float sumfvhelper2(v8sf x)
{
  v8sf t = __builtin_shuffle(x, (v8u32){4,5,6,7,0,1,2,3});
  x += t;
  t = __builtin_shuffle(x, (v8u32){2,3,0,1,6,7,4,5});
  x += t;
  t = __builtin_shuffle(x, (v8u32){1,0,3,2,5,4,7,6});
  x += t;
  return x[0];
}

float sumfv(float *x)
{
  //x = __builtin_assume_aligned(x, 64);
  v8sf *vx = (v8sf*)x;
  v8sf sumvv[4] = {{0}};
  for (size_t i = 0; i < 2048/8; i+=4)
    {
      sumvv[0] += vx[i+0];
      sumvv[1] += vx[i+1];
      sumvv[2] += vx[i+2];
      sumvv[3] += vx[i+3];
    }
  v8sf sumv = sumfvhelper1(sumvv);
  return sumfvhelper2(sumv);
}

哪个 gcc 4.8.4 gcc -Wall -Wextra -Wpedantic -std=gnu11 -march=native -O3 -fno-signed-zeros -fno-trapping-math -freciprocal-math -ffinite-math-only -fassociative-math -S 变成：

sumfv:
    vxorps  %xmm2, %xmm2, %xmm2
    xorl    %eax, %eax
    vmovaps %ymm2, %ymm3
    vmovaps %ymm2, %ymm0
    vmovaps %ymm2, %ymm1
.L7:
    addq    $4, %rax
    vaddps  (%rdi), %ymm1, %ymm1
    subq    $-128, %rdi
    vaddps  -96(%rdi), %ymm0, %ymm0
    vaddps  -64(%rdi), %ymm3, %ymm3
    vaddps  -32(%rdi), %ymm2, %ymm2
    cmpq    $256, %rax
    jne .L7
    vaddps  %ymm2, %ymm3, %ymm2
    vaddps  %ymm0, %ymm1, %ymm0
    vaddps  %ymm0, %ymm2, %ymm0
    vperm2f128  $1, %ymm0, %ymm0, %ymm1
    vaddps  %ymm0, %ymm1, %ymm0
    vpermilps   $78, %ymm0, %ymm1
    vaddps  %ymm0, %ymm1, %ymm0
    vpermilps   $177, %ymm0, %ymm1
    vaddps  %ymm0, %ymm1, %ymm0
    vzeroupper
    ret

第二个辅助函数并不是绝对必要的，但是对向量的元素求和往往会在 gcc 中产生糟糕的代码。如果你愿意做平台相关的内在函数，你可以用__builtin_ia32_hadps256()替换大部分。

【讨论】：