Comment activer sse3 autovectorization dans gcc

Question

J'ai une simple boucle avec le produit de n nombres complexes. Comme je réalise cette boucle des millions de fois, je veux que ce soit le plus rapide possible. Je comprends qu'il est possible de le faire rapidement en utilisant SSE3 et gcc intrinsics mais je suis intéressé à savoir s'il est possible d'obtenir gcc pour auto-vectoriser le code. Voici quelques exemples de code

#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; 
} 

L'ensemble vous obtenez de gcc -S-math -O3 -ffast est:

 .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

Le problème est que le type complex n'est pas compatible SIMD. Je n'ai jamais été un fan du type complex parce que c'est un objet composite qui ne correspond généralement pas à un type primitif ou une opération unique dans le matériel (certainement pas avec le matériel x86).

Afin de simplifier l'arithmétique SIMD, vous devez utiliser simultanément plusieurs numéros complexes. Pour SSE, vous devez opérer sur quatre numéros complexes à la fois.

Nous pouvons utiliser les codes vector extensions de GCC pour faciliter la syntaxe.

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

Ensuite, nous pouvons delcare une union d'un tableau et l'extension vecteur

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

Enfin, nous définissons un bloc de quatre nombres complexes comme celui-ci

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

où x est quatre parties réelles et y est quatre composants imaginaires.

Une fois que nous avons ce que nous pouvons plusieurs 4 nombres complexes à la fois comme celui-ci

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}; 
} 

et enfin nous arrivons à votre fonction modifiée pour fonctionner sur quatre nombres complexes à la fois.

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; 
} 

Regardons l'ensemble (syntaxe Intel) pour voir si elle est optimale

.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 

C'est exactement quatre 4 à l'échelle multiplications, un 4 à l'échelle plus, et une soustraction 4 à l'échelle. La variable p reste dans le registre et seul le tableau x est chargé de la mémoire comme nous le voulons.

Regardons à l'algèbre pour le produit des nombres complexes

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

C'est exactement quatre multiplications, une addition et une soustraction.

Comme je l'ai expliqué dans ce answer SIMD efficace algébriquement est souvent identique à l'arithmétique scalaire. Nous avons donc remplacé quatre multiplications de 1 largeur, addition et soustraction, avec quatre multiplications de 4, addition et soustraction. C'est le meilleur que vous pouvez faire avec SIMD 4-wide: quatre pour le prix d'un. Notez que cela n'a pas besoin d'instructions au-delà de SSE et aucune instruction SSE supplémentaire (à l'exception de FMA4) ne sera meilleure. Donc, sur un système 64 bits, vous pouvez compiler avec -O3.

Il est trivial d'étendre ceci pour SIMD 8-wide avec AVX.

Un avantage majeur de l'utilisation des extensions vectorielles de GCC est que vous obtenez FMA sans effort supplémentaire. Par exemple. si vous compilez avec -O3 -mfma4 la boucle principale est

.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 

Answer 2

Je ne suis pas un expert de montage mais j'ai réussi suivante. J'aurais des commentaires, mais il est trop grand:

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 

Ma commande de compilation était gcc -S -O3 -ffast-math -ftree-vectorizer-verbose=3 -ftree-slp-vectorize -ftree-vectorize -msse3 test.c vous n'avez pas besoin tous aussi peu obtient permis à -O3. Référez-vous à https://gcc.gnu.org/projects/tree-ssa/vectorization.html

Bien que je n'ai pas de réponse, j'ai essayé d'aider. Quand je mon architecture précise cpu (build) et je me suivant:

.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 

La commande est maintenant gcc -S -O3 -ffast-math -msse4 -march=haswell test.c où Haswell est mon cpu 4770HQ Core i7. Référez this pour votre cpu. Donc, comme vous voyez le jeu d'instructions AVX venir en image dans la deuxième version.

Un point de référence de l'échantillon de code ci-dessous:

$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; 
} 

Le réseau est statique si la majeure partie est-à-dire sur le disque ssd disque dur. Vous pouvez l'allouer en mémoire pour un traitement encore plus rapide. C'est des boucles de 200M. Binaire est 1.5G donc la plupart du temps est IO. Le CPU l'éclaire même sans -msse3 et -march. Tout ce dont vous avez besoin est -fast-math. Cela provoque une grande différence.

J'ai changé le programme suivant:

#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; 
} 

et compilé avec gcc -S -O3 -ffast-math -msse3 -mfpmath=sse -mavx -march=haswell test.c qui se traduit par:

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 

Alors ce qui semble me est que pour forcer gcc utiliser SSE3 vous noeud au code d'une certaine manière. http://sci.tuomastonteri.fi/programming/sse vous sera utile. Remarques finales: Si vous expérimentez avec différentes valeurs de limite supérieure pour i, vous verrez que différentes instructions sont produites. Je pense que la raison en est que gcc n'évalue pas les variables, donc vous pouvez utiliser des templates C++ capables de calculer le temps de compilation et de le faire.