如何在 gcc 中为复数启用 SSE3 addsubps 自动向量化？

Question

我有一个简单的循环，取 n 个复数的乘积。当我执行这个循环数百万次时，我希望它尽可能快。我知道可以使用 SSE3 和 gcc 内在函数快速完成此操作_mm_addsub_ps，但我感兴趣的是是否可以让 gcc 自动向量化这样的代码，即复数的乘积：

#include <complex.h>
complex float f(complex float x[], int n ) {
  complex float p = 1.0;
  for (int i = 0; i < n; i++)
    p *= x[i];
  return p;
}

您获得的程序集gcc -S -O3 -ffast-math是：

        .file   "test.c"
        .section        .text.unlikely,"ax",@progbits
.LCOLDB2:
        .text
.LHOTB2:
        .p2align 4,,15
        .globl  f
        .type   f, @function
f:
.LFB0:
        .cfi_startproc
        testl   %esi, %esi
        jle     .L4
        leal    -1(%rsi), %eax
        pxor    %xmm2, %xmm2
        movss   .LC1(%rip), %xmm3
        leaq    8(%rdi,%rax,8), %rax
        .p2align 4,,10
        .p2align 3
.L3:
        movaps  %xmm3, %xmm5
        movaps  %xmm3, %xmm4
        movss   (%rdi), %xmm0
        addq    $8, %rdi
        movss   -4(%rdi), %xmm1
        mulss   %xmm0, %xmm5
        mulss   %xmm1, %xmm4
        cmpq    %rdi, %rax
        mulss   %xmm2, %xmm0
        mulss   %xmm2, %xmm1
        movaps  %xmm5, %xmm3
        movaps  %xmm4, %xmm2
        subss   %xmm1, %xmm3
        addss   %xmm0, %xmm2
        jne     .L3
        movaps  %xmm2, %xmm1
.L2:
        movss   %xmm3, -8(%rsp)
        movss   %xmm1, -4(%rsp)
        movq    -8(%rsp), %xmm0
        ret
.L4:
        movss   .LC1(%rip), %xmm3
        pxor    %xmm1, %xmm1
        jmp     .L2
        .cfi_endproc
.LFE0:
        .size   f, .-f
        .section        .text.unlikely
.LCOLDE2:
        .text
.LHOTE2:
        .section        .rodata.cst4,"aM",@progbits,4
        .align 4
.LC1:
        .long   1065353216
        .ident  "GCC: (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609"
        .section        .note.GNU-stack,"",@progbits

Answer 1

问题在于该complex类型对 SIMD 不友好。我从来都不是该类型的粉丝complex，因为它是一个复合对象，通常不会映射到硬件中的原始类型或单个操作（当然不是 x86 硬件）。

为了使复杂算术 SIMD 友好，您需要同时操作多个复数。对于 SSE，您需要同时操作四个复数。

我们可以使用 GCC 的向量扩展来使语法更简单。

typedef float v4sf __attribute__ ((vector_size (16)));

然后我们可以删除数组和向量扩展的并集

typedef union {
  v4sf v;
  float e[4];
} float4

最后我们定义一个由四个复数组成的块，如下所示

typedef struct {
  float4 x;
  float4 y;
} complex4;

其中x是四个实部，y是四个虚部。

一旦我们有了这个，我们就可以像这样同时乘以 4 个复数

static complex4 complex4_mul(complex4 a, complex4 b) {
  return (complex4){a.x.v*b.x.v -a.y.v*b.y.v, a.y.v*b.x.v + a.x.v*b.y.v};
}

最后我们将您的函数修改为一次对四个复数进行操作。

complex4 f4(complex4 x[], int n) {
  v4sf one = {1,1,1,1};
  complex4 p = {one,one};
  for (int i = 0; i < n; i++) p = complex4_mul(p, x[i]);
  return p;
}

让我们看一下汇编（Intel 语法），看看它是否是最佳的

.L3:
    movaps  xmm4, XMMWORD PTR [rsi]
    add     rsi, 32
    movaps  xmm1, XMMWORD PTR -16[rsi]
    cmp     rdx, rsi
    movaps  xmm2, xmm4
    movaps  xmm5, xmm1
    mulps   xmm1, xmm3
    mulps   xmm2, xmm3
    mulps   xmm5, xmm0
    mulps   xmm0, xmm4
    subps   xmm2, xmm5
    addps   xmm0, xmm1
    movaps  xmm3, xmm2
    jne     .L3

这正好是四次 4 宽乘法、一次 4 宽加法和一次 4 宽减法。变量p保留在寄存器中，只有数组x从内存中加载，就像我们想要的那样。

让我们看一下复数乘积的代数

{a, bi}*{c, di} = {(ac - bd),(bc + ad)i}

这正好是四次乘法、一次加法和一次减法。

正如我在这个答案中所解释的，高效的 SIMD 在代数上通常与标量算术相同。因此，我们用四个 4 宽乘法、加法和减法替换了四个 1 宽乘法、加法和减法。这是使用 4 宽 SIMD 所能做到的最好效果：用 1 个的价格购买 4 个 SIMD。

请注意，这不需要 SSE 之外的任何指令，并且没有额外的 SSE 指令（FMA4 除外）会更好。因此在 64 位系统上您可以使用-O3.

使用 AVX 将其扩展为 8 宽 SIMD 是微不足道的。

使用 GCC 向量扩展的一大优势是您无需任何额外努力即可获得 FMA。例如，如果您使用-O3 -mfma4主循环进行编译是

.L3:
    vmovaps xmm0, XMMWORD PTR 16[rsi]
    add     rsi, 32
    vmulps  xmm1, xmm0, xmm2
    vmulps  xmm0, xmm0, xmm3
    vfmsubps        xmm1, xmm3, XMMWORD PTR -32[rsi], xmm1
    vmovaps xmm3, xmm1
    vfmaddps        xmm2, xmm2, XMMWORD PTR -32[rsi], xmm0
    cmp     rdx, rsi
    jne     .L3