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

qvk1mo1f  于 2023-01-09  发布在  其他
关注(0)|答案(2)|浏览(121)

我有一个简单的循环,需要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
a5g8bdjr

a5g8bdjr1#

问题是complex类型不是SIMD友好的,我从来都不是complex类型的粉丝,因为它是一个复合对象,通常不会Map到硬件中的原语类型或单个操作(当然在x86硬件中不会)。
为了使复数运算对SIMD友好,你需要同时操作多个复数。对于SSE,你需要同时操作四个复数。
我们可以使用GCC的vector extensions来简化语法。

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进制乘法,1个4进制加法和1个4进制减法,变量p保留在寄存器中,只有数组x从内存中加载,就像我们想要的那样。
让我们看看复数乘积的代数

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

这就是四次乘法,一次加法,一次减法。
正如我在本答案中所解释的,高效的SIMD代数运算通常与标量算术运算相同。因此,我们将4个1宽度的乘法、加法和减法替换为4个4宽度的乘法、加法和减法。这是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
lsmd5eda

lsmd5eda2#

我不是一个装配Maven,但我已经管理以下。我会评论,但它太大:

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,其中haswell是我的i74770hq cpu。请参考this了解你的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。CPU即使没有-msse3和-march也会很忙。你所需要的只是-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++模板,它能够计算编译时间。

相关问题