如何要求GCC完全展开这个循环（即剥离这个循环）？

Question

是否有一种方法可以指示GCC（我使用的是4.8.4）完全展开底部函数中的＜code＞循环，即剥离该循环？循环的迭代次数在编译时是已知的：58。

让我先解释一下我所做的努力。

通过检查GAS输出：

gcc -fpic -O2 -S GEPDOT.c

使用了12个寄存器XMM0 - XMM11。如果我将标志< code>-funroll-loops传递给gcc:

gcc -fpic -O2 -funroll-loops -S GEPDOT.c

循环只展开两次。我检查了GCC优化选项。GCC说-funrol-loops也会打开-frename-ynster，所以当GCC展开循环时，它在寄存器分配方面的首选是使用“剩余”寄存器。但是XMM12-XMM15只剩下4个，所以GCC最好只能展开2次。如果有48个而不是16个XMM寄存器可用，GCC将毫不费力地展开4次。

然而我做了另一个实验。我首先手动展开while循环两次，获得一个函数< code>GEPDOT_2。那么两者之间就没有任何区别

gcc -fpic -O2 -S GEPDOT_2.c

和

gcc -fpic -O2 -funroll-loops -S GEPDOT_2.c

由于＜code＞GEPDOT_2＜/code＞已用完所有寄存器，因此不执行展开。

GCC确实注册重命名以避免引入潜在的错误依赖。但我肯定知道在我的GEPDOT中不会有这样的潜在性；即使有，也不重要。我自己尝试了展开循环，展开4次比展开2次快，比不展开快。当然我可以手动展开更多次，但很繁琐。GCC能帮我做这个吗？谢谢。

// C file "GEPDOT.c"
#include <emmintrin.h>

void GEPDOT (double *A, double *B, double *C) {
  __m128d A1_vec = _mm_load_pd(A); A += 2;
  __m128d B_vec = _mm_load1_pd(B); B++;
  __m128d C1_vec = A1_vec * B_vec;
  __m128d A2_vec = _mm_load_pd(A); A += 2;
  __m128d C2_vec = A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  __m128d C3_vec = A1_vec * B_vec;
  __m128d C4_vec = A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  __m128d C5_vec = A1_vec * B_vec;
  __m128d C6_vec = A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  __m128d C7_vec = A1_vec * B_vec;
  A1_vec = _mm_load_pd(A); A += 2;
  __m128d C8_vec = A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  int k = 58;
  /* can compiler unroll the loop completely (i.e., peel this loop)? */
  while (k--) {
    C1_vec += A1_vec * B_vec;
    A2_vec = _mm_load_pd(A); A += 2;
    C2_vec += A2_vec * B_vec;
    B_vec = _mm_load1_pd(B); B++;
    C3_vec += A1_vec * B_vec;
    C4_vec += A2_vec * B_vec;
    B_vec = _mm_load1_pd(B); B++;
    C5_vec += A1_vec * B_vec;
    C6_vec += A2_vec * B_vec;
    B_vec = _mm_load1_pd(B); B++;
    C7_vec += A1_vec * B_vec;
    A1_vec = _mm_load_pd(A); A += 2;
    C8_vec += A2_vec * B_vec;
    B_vec = _mm_load1_pd(B); B++;
    }
  C1_vec += A1_vec * B_vec;
  A2_vec = _mm_load_pd(A);
  C2_vec += A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  C3_vec += A1_vec * B_vec;
  C4_vec += A2_vec * B_vec;
  B_vec = _mm_load1_pd(B); B++;
  C5_vec += A1_vec * B_vec;
  C6_vec += A2_vec * B_vec;
  B_vec = _mm_load1_pd(B);
  C7_vec += A1_vec * B_vec;
  C8_vec += A2_vec * B_vec;
  /* [write-back] */
  A1_vec = _mm_load_pd(C); C1_vec = A1_vec - C1_vec;
  A2_vec = _mm_load_pd(C + 2); C2_vec = A2_vec - C2_vec;
  A1_vec = _mm_load_pd(C + 4); C3_vec = A1_vec - C3_vec;
  A2_vec = _mm_load_pd(C + 6); C4_vec = A2_vec - C4_vec;
  A1_vec = _mm_load_pd(C + 8); C5_vec = A1_vec - C5_vec;
  A2_vec = _mm_load_pd(C + 10); C6_vec = A2_vec - C6_vec;
  A1_vec = _mm_load_pd(C + 12); C7_vec = A1_vec - C7_vec;
  A2_vec = _mm_load_pd(C + 14); C8_vec = A2_vec - C8_vec;
  _mm_store_pd(C,C1_vec); _mm_store_pd(C + 2,C2_vec);
  _mm_store_pd(C + 4,C3_vec); _mm_store_pd(C + 6,C4_vec);
  _mm_store_pd(C + 8,C5_vec); _mm_store_pd(C + 10,C6_vec);
  _mm_store_pd(C + 12,C7_vec); _mm_store_pd(C + 14,C8_vec);
  }

感谢@user3386109的评论，我想对这个问题进行一点扩展@user3386109提出了一个非常好的问题。实际上，当有这么多并行指令需要调度时，我确实对编译器的最佳寄存器分配能力有一些怀疑。

我个人认为一个可靠的方法是首先在asm内联汇编中编写循环体(这是HPC的关键),然后根据我的需要多次复制它。今年早些时候，我发表了一篇不受欢迎的文章:内嵌汇编。代码略有不同，因为循环迭代次数j是一个函数参数，因此在编译时是未知的。在这种情况下，我不能完全展开循环，所以我只复制了两次汇编代码，并将循环转换为标签和跳转。结果是，我编写的汇编程序的性能比编译器生成的汇编程序高5%,这可能意味着编译器无法以我们期望的最佳方式分配寄存器。

我是（现在仍然是）汇编编程的一个婴儿，所以这是一个很好的案例研究，可以让我学习一点x86汇编。但从长远来看，我不倾向于代码＜code＞GEPDOT

1. asm内联程序集因不可移植而受到批评。虽然我不明白为什么。也许是因为不同的机器有不同的寄存器被破坏
2. 编译器也越来越好。所以，我仍然更喜欢算法优化和更好的C编码习惯，以帮助编译器生成良好的输出
3. 最后一个原因更重要。迭代次数可能不总是58。我正在开发一个高性能矩阵分解子程序。对于缓存阻塞因子nb，迭代次数将为nb-2。我不会把＜code＞nb

如果你也能分享一些生产高性能、可移植库的经验，我将非常感激。

Answer 1

尝试调整优化器参数：

gcc -O3 -funroll-loops --param max-completely-peeled-insns=1000 --param max-completely-peel-times=100

这应该能行。

Answer 2

这不是一个答案，但其他试图用GCC对矩阵乘法进行矢量化的人可能会感兴趣。

下面，我假设c是行主顺序的4×4矩阵，a是列主顺序的4行n列矩阵(转置)，b是行主顺序的4列n行矩阵，要计算的运算是c = a × b c，其中×表示矩阵乘法。

完成此操作的幼稚功能是

void slow_4(double       *c,
            const double *a,
            const double *b,
            size_t        n)
{
    size_t row, col, i;

    for (row = 0; row < 4; row++)
        for (col = 0; col < 4; col++)
            for (i = 0; i < n; i++)
                c[4*row+col] += a[4*i+row] * b[4*i+col];
}

GCC使用为SSE2 / SSE3生成了非常好的代码

#if defined(__SSE2__) || defined(__SSE3__)

typedef  double  vec2d  __attribute__((vector_size (2 * sizeof (double))));

void fast_4(vec2d       *c,
            const vec2d *a,
            const vec2d *b,
            size_t       n)
{
    const vec2d *const b_end = b + 2L * n;

    vec2d s00 = c[0];
    vec2d s02 = c[1];
    vec2d s10 = c[2];
    vec2d s12 = c[3];
    vec2d s20 = c[4];
    vec2d s22 = c[5];
    vec2d s30 = c[6];
    vec2d s32 = c[7];

    while (b < b_end) {
        const vec2d b0 = b[0];
        const vec2d b2 = b[1];
        const vec2d a0 = { a[0][0], a[0][0] };
        const vec2d a1 = { a[0][1], a[0][1] };
        const vec2d a2 = { a[1][0], a[1][0] };
        const vec2d a3 = { a[1][1], a[1][1] };
        s00 += a0 * b0;
        s10 += a1 * b0;
        s20 += a2 * b0;
        s30 += a3 * b0;
        s02 += a0 * b2;
        s12 += a1 * b2;
        s22 += a2 * b2;
        s32 += a3 * b2;
        b += 2;
        a += 2;
    }

    c[0] = s00;
    c[1] = s02;
    c[2] = s10;
    c[3] = s12;
    c[4] = s20;
    c[5] = s22;
    c[6] = s30;
    c[7] = s32; 
}

#endif

对于AVX，GCC可以做得更好

#if defined(__AVX__) || defined(__AVX2__)

typedef  double  vec4d  __attribute__((vector_size (4 * sizeof (double))));

void fast_4(vec4d       *c,
            const vec4d *a,
            const vec4d *b,
            size_t       n)
{
    const vec4d *const b_end = b + n;

    vec4d s0 = c[0];
    vec4d s1 = c[1];
    vec4d s2 = c[2];
    vec4d s3 = c[3];

    while (b < b_end) {
        const vec4d bc = *(b++);
        const vec4d ac = *(a++);
        const vec4d a0 = { ac[0], ac[0], ac[0], ac[0] };
        const vec4d a1 = { ac[1], ac[1], ac[1], ac[1] };
        const vec4d a2 = { ac[2], ac[2], ac[2], ac[2] };
        const vec4d a3 = { ac[3], ac[3], ac[3], ac[3] };
        s0 += a0 * bc;
        s1 += a1 * bc;
        s2 += a2 * bc;
        s3 += a3 * bc;
    }

    c[0] = s0;
    c[1] = s1;
    c[2] = s2;
    c[3] = s3;
}

#endif

使用gcc-4 . 8 . 4(< code >-O2-March = x86-64-mtune = generic-msse 3 )生成的程序集的SSE3版本本质上是

fast_4:
        salq    $5, %rcx
        movapd  (%rdi), %xmm13
        addq    %rdx, %rcx
        cmpq    %rcx, %rdx
        movapd  16(%rdi), %xmm12
        movapd  32(%rdi), %xmm11
        movapd  48(%rdi), %xmm10
        movapd  64(%rdi), %xmm9
        movapd  80(%rdi), %xmm8
        movapd  96(%rdi), %xmm7
        movapd  112(%rdi), %xmm6
        jnb     .L2
.L3:
        movddup (%rsi), %xmm5
        addq    $32, %rdx
        movapd  -32(%rdx), %xmm1
        addq    $32, %rsi
        movddup -24(%rsi), %xmm4
        movapd  %xmm5, %xmm14
        movddup -16(%rsi), %xmm3
        movddup -8(%rsi), %xmm2
        mulpd   %xmm1, %xmm14
        movapd  -16(%rdx), %xmm0
        cmpq    %rdx, %rcx
        mulpd   %xmm0, %xmm5
        addpd   %xmm14, %xmm13
        movapd  %xmm4, %xmm14
        mulpd   %xmm0, %xmm4
        addpd   %xmm5, %xmm12
        mulpd   %xmm1, %xmm14
        addpd   %xmm4, %xmm10
        addpd   %xmm14, %xmm11
        movapd  %xmm3, %xmm14
        mulpd   %xmm0, %xmm3
        mulpd   %xmm1, %xmm14
        mulpd   %xmm2, %xmm0
        addpd   %xmm3, %xmm8
        mulpd   %xmm2, %xmm1
        addpd   %xmm14, %xmm9
        addpd   %xmm0, %xmm6
        addpd   %xmm1, %xmm7
        ja      .L3
.L2:
        movapd  %xmm13, (%rdi)
        movapd  %xmm12, 16(%rdi)
        movapd  %xmm11, 32(%rdi)
        movapd  %xmm10, 48(%rdi)
        movapd  %xmm9, 64(%rdi)
        movapd  %xmm8, 80(%rdi)
        movapd  %xmm7, 96(%rdi)
        movapd  %xmm6, 112(%rdi)
        ret

生成的程序集的 AVX 版本 （-O2 -march=x86-64 -mtune=generic -mavx） 本质上是

fast_4:
        salq       $5, %rcx
        vmovapd    (%rdi), %ymm5
        addq       %rdx, %rcx
        vmovapd    32(%rdi), %ymm4
        cmpq       %rcx, %rdx
        vmovapd    64(%rdi), %ymm3
        vmovapd    96(%rdi), %ymm2
        jnb        .L2
.L3:
        addq       $32, %rsi
        vmovapd    -32(%rsi), %ymm1
        addq       $32, %rdx
        vmovapd    -32(%rdx), %ymm0
        cmpq       %rdx, %rcx
        vpermilpd  $0, %ymm1, %ymm6
        vperm2f128 $0, %ymm6, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm6, %ymm6
        vaddpd     %ymm6, %ymm5, %ymm5
        vpermilpd  $15, %ymm1, %ymm6
        vperm2f128 $0, %ymm6, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm6, %ymm6
        vaddpd     %ymm6, %ymm4, %ymm4
        vpermilpd  $0, %ymm1, %ymm6
        vpermilpd  $15, %ymm1, %ymm1
        vperm2f128 $17, %ymm6, %ymm6, %ymm6
        vperm2f128 $17, %ymm1, %ymm1, %ymm1
        vmulpd     %ymm0, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm1, %ymm0
        vaddpd     %ymm6, %ymm3, %ymm3
        vaddpd     %ymm0, %ymm2, %ymm2
        ja         .L3
.L2:
        vmovapd    %ymm5, (%rdi)
        vmovapd    %ymm4, 32(%rdi)
        vmovapd    %ymm3, 64(%rdi)
        vmovapd    %ymm2, 96(%rdi)
        vzeroupper
        ret

我猜寄存器调度不是最佳的，但看起来也不差。我个人对以上内容很满意，没有试图在这一点上手动优化它。

在 Core i5-4200U 处理器（支持 AVX2）上，上述函数的快速版本计算 1843 个 CPU 周期（中位数）和 AVX2 的 1248 个周期中两个 4×256 矩阵的乘积。这归结为每个矩阵条目 1.8 和 1.22 个周期。用于比较，未矢量化的慢速版本每个矩阵条目大约需要 11 个周期。

(循环计数是中间值，即一半的测试更快。我只运行了大约100k次重复的粗略基准测试，所以不要完全相信这些数字。)

在这个CPU上，缓存效果是这样的，4×512矩阵大小的AVX2仍然是每个条目1.2个周期，但在4×1024时，它下降到1.4,4×4096到1.5,4×8192到1.8，每个条目4×65536到2.2个周期。SSE3版本保持每个条目1.8个周期，直到4×3072，此时它开始变慢；在4×65536时，它也是每个条目大约2.2个周期。我确实相信这个（笔记本电脑！）CPU在这一点上是缓存带宽有限的。