4

我有一个简单的循环，取 n 个复数的乘积。当我执行这个循环数百万次时，我希望它尽可能快。我知道使用 SSE3 和 gcc 内在函数可以快速做到这一点，但我对是否可以让 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
4

2 回答 2

2

问题是该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 可以做到的最好的事情：一个价格为四个。

请注意，这不需要 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
于 2017-01-03T09:43:33.517 回答
1

我不是装配专家，但我已经成功了。我会评论但它太大了：

cat test.s
    .file   "test.c"
    .text
    .p2align 4,,15
    .globl  f
    .type   f, @function
f:
.LFB0:
    .cfi_startproc
    testl   %esi, %esi
    jle     .L4
    leal    -1(%rsi), %eax
    pxor    %xmm0, %xmm0
    movss   .LC1(%rip), %xmm1
    leaq    8(%rdi,%rax,8), %rax
    .p2align 4,,10
    .p2align 3
.L3:
    movaps  %xmm1, %xmm4
    movss   (%rdi), %xmm3
    movss   4(%rdi), %xmm2
    mulss   %xmm3, %xmm1
    mulss   %xmm2, %xmm4
    addq    $8, %rdi
    mulss   %xmm0, %xmm2
    cmpq    %rdi, %rax
    mulss   %xmm3, %xmm0
    subss   %xmm2, %xmm1
    addss   %xmm4, %xmm0
    jne     .L3
.L1:
    movss   %xmm1, -8(%rsp)
    movss   %xmm0, -4(%rsp)
    movq    -8(%rsp), %xmm0
    ret
.L4:
    movss   .LC1(%rip), %xmm1
    pxor    %xmm0, %xmm0
    jmp     .L1
    .cfi_endproc
.LFE0:
    .size   f, .-f
    .section        .rodata.cst4,"aM",@progbits,4
    .align 4
.LC1:
    .long   1065353216
    .ident  "GCC: (Ubuntu 6.2.0-5ubuntu12) 6.2.0 20161005"
    .section        .note.GNU-stack,"",@progbits

我的编译命令是gcc -S -O3 -ffast-math -ftree-vectorizer-verbose=3 -ftree-slp-vectorize -ftree-vectorize -msse3 test.c您不需要所有这些，因为在 -O3 启用的很少。参考https://gcc.gnu.org/projects/tree-ssa/vectorization.html

虽然我没有答案，但我试图提供帮助。当我指定我的 cpu 架构（构建）时，我得到以下信息：

    .file   "test.c"
    .text
    .p2align 4,,15
    .globl  f
    .type   f, @function
f:
.LFB0:
    .cfi_startproc
    testl   %esi, %esi
    jle     .L4
    vmovss  .LC1(%rip), %xmm1
    leal    -1(%rsi), %eax
    vxorps  %xmm0, %xmm0, %xmm0
    leaq    8(%rdi,%rax,8), %rax
    .p2align 4,,10
    .p2align 3
.L3:
    vmovss  (%rdi), %xmm2
    vmovss  4(%rdi), %xmm3
    addq    $8, %rdi
    vmulss  %xmm3, %xmm0, %xmm4
    vmulss  %xmm2, %xmm0, %xmm0
    vfmadd231ss     %xmm3, %xmm1, %xmm0
    vfmsub132ss     %xmm2, %xmm4, %xmm1
    cmpq    %rdi, %rax
    jne     .L3
.L1:
    vmovss  %xmm1, -8(%rsp)
    vmovss  %xmm0, -4(%rsp)
    vmovq   -8(%rsp), %xmm0
    ret
.L4:
    vmovss  .LC1(%rip), %xmm1
    vxorps  %xmm0, %xmm0, %xmm0
    jmp     .L1
    .cfi_endproc
.LFE0:
    .size   f, .-f
    .section        .rodata.cst4,"aM",@progbits,4
    .align 4
.LC1:
    .long   1065353216
    .ident  "GCC: (Ubuntu 6.2.0-5ubuntu12) 6.2.0 20161005"
    .section        .note.GNU-stack,"",@progbits

现在的命令gcc -S -O3 -ffast-math -msse4 -march=haswell test.c是我的 i7 4770HQ cpu 所在的 haswell。请参阅此为您的 cpu。

因此，正如您在第二个版本中看到的 AVX 指令集一样。

以下代码的示例基准：

$time ./a.out 
0.000000
real    0m0.684s
user    0m0.620s
sys     0m0.060s


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

int main()
{
    static complex float x[200000000] = {0.0, 1.0, 2.0, 4.0, 5.0, 6.0};
    complex float p = f(x, 200000000);

    printf("%f", creal(p));

    return 0;
}

该阵列是静态的，因此大部分都在磁盘上，即 ssd 硬盘驱动器上。您可以将其分配到内存中以加快处理速度。这是 200M 循环。二进制是1.5G，所以大部分时间都是IO。即使没有 -msse3 和 -march，CPU 也很出色。您只需要-ffast-math。这造成了很大的不同。

我将程序更改为以下内容：

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

int main() {
    float x[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};

    printf("%f\n", f(x, 8));

    return 0;
}

并编译gcc -S -O3 -ffast-math -msse3 -mfpmath=sse -mavx -march=haswell test.c结果：

f:
.LFB23:
    .cfi_startproc
    vmovups (%rdi), %ymm2
    vxorps  %xmm1, %xmm1, %xmm1
    vperm2f128      $33, %ymm1, %ymm2, %ymm0
    vmulps  %ymm2, %ymm0, %ymm0
    vperm2f128      $33, %ymm1, %ymm0, %ymm2
    vshufps $78, %ymm2, %ymm0, %ymm2
    vmulps  %ymm2, %ymm0, %ymm0
    vperm2f128      $33, %ymm1, %ymm0, %ymm1
    vpalignr        $4, %ymm0, %ymm1, %ymm1
    vmulps  %ymm1, %ymm0, %ymm0
    vzeroupper
    ret
    .cfi_endproc

所以在我看来，要强制 gcc 使用 SSE3 你的节点以某种方式编码。http://sci.tuomastonteri.fi/programming/sse将对您有用。

最后说明：如果您尝试 i 的不同上限值，您将看到产生了不同的指令。我认为这样做的原因是 gcc 不评估变量，因此您可能希望使用能够编译时间计算并执行此操作的 C++ 模板。

于 2017-01-01T16:59:03.953 回答