如何在 gcc 中启用 sse3 自动向量化
How to enable sse3 autovectorization in gcc
我有一个简单的循环，它取 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 , %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
我不是装配专家，但我已经成功了。我会发表评论，但它太大了：
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 , %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 , %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 是我的 i7 4770HQ 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。即使没有 -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 , %ymm1, %ymm2, %ymm0
vmulps %ymm2, %ymm0, %ymm0
vperm2f128 , %ymm1, %ymm0, %ymm2
vshufps , %ymm2, %ymm0, %ymm2
vmulps %ymm2, %ymm0, %ymm0
vperm2f128 , %ymm1, %ymm0, %ymm1
vpalignr , %ymm0, %ymm1, %ymm1
vmulps %ymm1, %ymm0, %ymm0
vzeroupper
ret
.cfi_endproc
所以在我看来，要强制 gcc 使用 SSE3，您需要以某种方式节点编码。 http://sci.tuomastonteri.fi/programming/sse 对你有用。
最后说明：如果您对 i 的上限值使用不同的值进行试验，您会看到生成了不同的指令。我认为这样做的原因是 gcc 不计算变量，所以你可能想使用能够编译时计算的 C++ 模板并执行它。
问题是 complex 类型对 SIMD 不友好。我从来都不是 complex 类型的粉丝，因为它是一个复合对象，通常不会映射到原始类型或硬件中的单一操作（当然不是 x86 硬件）。
为了使复杂算术 SIMD 友好，您需要同时对多个复数进行运算。对于 SSE，您需要一次对四个复数进行运算。
我们可以使用 GCC 的 vector extensions 来简化语法。
typedef float v4sf __attribute__ ((vector_size (16)));
然后我们可以delcare一个数组和向量扩展的并集
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}
正好是四次乘法、一次加法和一次减法。
正如我在此 answer 中解释的那样，高效的 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
我有一个简单的循环，它取 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 , %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
我不是装配专家，但我已经成功了。我会发表评论，但它太大了：
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 , %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 , %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 是我的 i7 4770HQ 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。即使没有 -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 , %ymm1, %ymm2, %ymm0
vmulps %ymm2, %ymm0, %ymm0
vperm2f128 , %ymm1, %ymm0, %ymm2
vshufps , %ymm2, %ymm0, %ymm2
vmulps %ymm2, %ymm0, %ymm0
vperm2f128 , %ymm1, %ymm0, %ymm1
vpalignr , %ymm0, %ymm1, %ymm1
vmulps %ymm1, %ymm0, %ymm0
vzeroupper
ret
.cfi_endproc
所以在我看来，要强制 gcc 使用 SSE3，您需要以某种方式节点编码。 http://sci.tuomastonteri.fi/programming/sse 对你有用。
最后说明：如果您对 i 的上限值使用不同的值进行试验，您会看到生成了不同的指令。我认为这样做的原因是 gcc 不计算变量，所以你可能想使用能够编译时计算的 C++ 模板并执行它。
问题是 complex 类型对 SIMD 不友好。我从来都不是 complex 类型的粉丝，因为它是一个复合对象，通常不会映射到原始类型或硬件中的单一操作（当然不是 x86 硬件）。
为了使复杂算术 SIMD 友好，您需要同时对多个复数进行运算。对于 SSE，您需要一次对四个复数进行运算。
我们可以使用 GCC 的 vector extensions 来简化语法。
typedef float v4sf __attribute__ ((vector_size (16)));
然后我们可以delcare一个数组和向量扩展的并集
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}
正好是四次乘法、一次加法和一次减法。
正如我在此 answer 中解释的那样，高效的 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