MMX, SSE, AVX: tout sur la programmation SIMD, la vectorisation, les intrinsics

[vectra]

Ou comment vectra lance un topic de vectorisation


Bon alors, la vectorisation c'est le bien, alors il faut en mettre partout

Ce petit topic pour rassembler un minimum de documentation et de tutoriels sur la question, et aussi poser vos questions aux Men In Black sans pour autant dégoûter le reste des tauliers du topic de la programmation.


In a nutshell: la plupart des CPU Desktop possèdent dorénavant des unités de calcul vectorisées, qui ont remplacé presque toutes les unités scalaires depuis genre 1999 (pour le SSE). Maintenant, le scalaire est implémenté par une sous-partie du vectorisé, donc autant y aller franco.


Alors les unités vectorisées, qu'est-ce à dire? Visiblement, les fondeurs n'ont rien trouvé de plus drôle à faire que de rajouter des séries d'additionneurs, de multiplicateurs (etc) dans les unités d'addition, de multiplication (etc). Et donc, cette chaine d'unités peut fonctionner en simultané, chaque sous-unité ne traitant que les deux opérandes dont elle a la charge.

En conséquence, vous pouvez programmer une série d'opérations similaires qui vont s'effectuer en une seule passe sur des séries de nombres. L'addition flottante, par exemple, prend 3 cycles CPU pour traiter deux nombres flottants: on parle alors d'opération scalaires et d'opérandes scalaires, soit le mode par défaut quand vous écrivez un programme.
Hé bien figurez vous que l'addition vectorielle, traitant par exemple deux vecteurs de 4 flottants chacun, ne va prendre que 3 cycles CPU. Pas 12: 3.


Les additions se font en parallèle. Donc ouais, potentiellement, on peut réduire la durée en cycles CPU des calculs arithmétiques (entre autres) d'un facteur égal à la longueur du vecteur.

Bon, on a vu que les unités vectorisées, c'était du concret, du matos. Mais ces vecteurs alors? Hé bien, Intel et les autres rigolos ont fourni des jeux d'instructions assembleur et de registres assembleur qui permettent de programmer directement en vectorisé. On verra comment laisser ça aux masochistes et aux trve elite en se contentant de programmer ça en C via les Intrinsics.
Parmi ces jeux, on distingue le MMX (l'historique, celui du Pentiom 1, dont tout le monde se foutait à l'époque), le SSE qui rigole moins (v1 à 4 au moins), et le récent AVX.Les instructions: on vient d'en voir une, implémentée par une unité vectorisée flottante. Ben le vecteur, il est implémenté comme un long registre dont un set est fourni avec le jeu disponible sur notre CPU. Pour info, on ne perd jamais de fonctionnalité sur un processeur Intel/AMD: si vous avez l'AVX, vous avez du SSE4/3/2/1 et du MMX, donc au final ça vous fait du peuple, et en registres, et en instructions assembleur/intrinsics.

Sur le SSE4, vos registres ne font "que" 128 bits, et vos unités vectorisées ne peuvent pas traiter plus grand. Ca vous permet déjà de stocker 4 flottants 32 bits en simple précision, 2 flottants 64 bits, 8 variables 16 bits, et 16 char. Sur l'AVX (dispo dès le Core i5 2500k en gros), ces registres passent à 256 bits: ça rigole déjà moins dans le fond, hein
Si j'ai bien compris, depuis Sandy jusqu'à Skylake, on reste encore en 256 bits: seul le jeu d'instruction a été progressivement étendu de notre point de vue. Par contre, le Xéon Phi pour les grands (Intel Mic, Knight's Corner, etc) est le premier à implémenter des registres et instructions en 512 bits. Ca doit faire pile une ligne de cache mémoire, si j'ai bien compris.


Les instructions c'est bien, mais la mémoire?

Rien ne sert d'accélérer encore l'arithmétique si vous faites beaucoup d'accès mémoire. Le CPU en scalaire est déjà très limité par la mémoire pour des calculs simples sur de grands flux de données, en audio, en vidéo, en grosse 2-D, en 3-D, etc. On parle alors de streaming.

C'est pourquoi beaucoup de tutoriaux parlent par exemple de fractales de Mandelbrot, où qu'on fait plein de calculs en arithmétique complexe pour décider de la couleur d'un pixel, avant de l'écrire en mémoire. Dans ce cas de figure, les SSE/AVX font merveille, et parviennent à réduire le temps de calcul qui occulte le temps d'accès à la mémoire.

Mais le SSE a aussi été l'occasion de régler ce problème très chiant de bande passante mémoire pour des algos multimédia (MMX à l'origine). Le SSE fournit des instructions pour contourner carrément le cache, ce qui peut par exemple doubler les perfs en écriture (ouais, sans déconner) par rapport à un code trivial. Le SSE permet aussi de s'adresser au prefetcher (...)


Des tutoriels pour commencer:

Insérer ici une liste de tutos
Alab m'a envoyé des candidats, ça arrive.



Des docs de référence pour continuer:

une liste exhaustive et bien lisible des intrinsics selon Intel: https://software.intel.com/sites/lan...trinsicsGuide/

En cas d'erreur trop flagrante, je corrigerai!

[Sp1d3r]

Quid de la vectorisation dans un contexte avec des indirections inévitables ?

[vectra]

Ben j'ai eu le cas il y a peu: normalement, c'est complètement DTC, mais bizarrement les perfs n'étaient pas si mauvaises (mais j'ai pas terminé l'implémentation).
Je posterai ça si tu veux la terminer

[Møgluglu]

Quid de la vectorisation dans un contexte avec des indirections inévitables ?

Tu as l'instruction Gather dans AVX 2, et les instructions Gather et Scatter dans AVX-512. Tu passes un pointeur scalaire de base x et un vecteur d'offsets i et tu récupère le vecteur composé des x[i]. C'est une seule instruction, mais ça peut être plus ou moins long suivant l'emplacement des i. Mais si tu en as plusieurs qui tombent sur la même ligne de cache, ça se passe plutôt bien.

D'ailleurs sur GPU on ne se pose même pas la question, on n'a que des Gather et Scatter, pas de Load / Store avec adresse scalaire.

[gregounech]

/me s'est abonné à ce topic.

Sinon, c'est quoi cette histoire d'Intrinsics ?

Des fonctions C qui indiquent au compilo que tout çe qui est fait dedans est vectorisable, genre tu te fais une fonction qui prends 2 tableaux de 4 floats et qui les ajoutent comme dans les exemples ?

[vectra]

Désolé pour le lag, je dois rendre un bout de manuscrit viteuf.
C'est con, parce qu'il y a justement une section SSE pour les nuls encadrants dedans.
Non, c'est carrément une périphrase pour manipuler directement une instruction assembleur, mais dans du code C, et sans avoir à compter ses registres x86.

[gregounech]

Non, c'est carrément une périphrase pour manipuler directement une instruction assembleur, mais dans du code C, et sans avoir à compter ses registres x86.

Ok, donc un espèce d'espace d'abstraction entre les instructions assembleurs et du code C classique. Fais moi cette instruction sur ces données, tu te démerdes pour les détails?

[vectra]

Non, c'est applique telle instruction assembleur sur ces opérandes, point.

C'est l'équivalent "facile à faire" de la programmation assembleur. C'est plutôt chiant de coder toute une fonction en assembleur, et c'est assez chiant aussi de faire des inclusions de code assembleur dans du C (en tous cas, avec gcc). Dans le second cas, tu dois associer des registres du langage x86 avec des variables avant et après, et ça n'a rien d'aisé je trouve.
D'autant plus que le langage x86 oblige à utiliser un petit nombre de registres de par la norme, alors que derrière, les CPU sont capables d'en gérer plusieurs centaines.
Perso, ce qui me tue, ce sont en plus les notations assembleur différentes entre Intel et les autres

Un intrinsic est une instruction C chopée au vol par le compilateur.
L'avantage d'utiliser les intrinsics pour vectoriser, c'est que t'es sûr que ton code sera bien et complètement vectorisé, vu qu'il sera écrit ainsi dès le départ. L'avantage de cette approche, c'est que tu peux coder en C standard, repérer les sections lentes dans le code, et éventuellement ne vectoriser à la main que la fonction ou le corps de boucle qui est sollicité.

Une grosse différence entre assembleur et intrinsics, c'est que dans le second cas, le compilateur (son Backend, si j'ai bien compris) est en mesure de réordonner et modifier le code x86 pour optimiser son fonctionnement (exemples wanted, ça va encore parler IACA et latences je sens). Alors que du code assembleur n'est jamais retouché, ce qui gâche des opportunités d'opti.

[gregounech]

Non, c'est applique telle instruction assembleur sur ces opérandes, point.

C'est l'équivalent "facile à faire" de la programmation assembleur. C'est plutôt chiant de coder toute une fonction en assembleur, et c'est assez chiant aussi de faire des inclusions de code assembleur dans du C (en tous cas, avec gcc). Dans le second cas, tu dois associer des registres du langage x86 avec des variables avant et après, et ça n'a rien d'aisé je trouve.
D'autant plus que le langage x86 oblige à utiliser un petit nombre de registres de par la norme, alors que derrière, les CPU sont capables d'en gérer plusieurs centaines.
Perso, ce qui me tue, ce sont en plus les notations assembleur différentes

Okay, j'ai jamais fait d'inclusions de code assembleur dans du C en x86, seulement en MIPS dans le cadre de cours en se servant de QEMU et dans des cas très limités / petits et avec un squelette déja très largement fourni.

[vectra]

Je me débrouillais en MIPS, mais en x86, je suis un peu une grosse quiche.

Un truc marrant: tu peux voir le code assembleur en tout joli sur cette appli en ligne:
https://gcc.godbolt.org/

C'est bien plus lisible que ce que j'arrive à faire cracher à gcc.

Typiquement, mulss opinout opin => multiplication normale deux scalaires flottants, et mulps vectorielle entre deux vecteurs (no way ).
Fun fact: même en scalaire, tu as un ou deux vecteurs comme opérandes, étant donné que cette instruction est traitée par la même unité vectorisée (e.g AVX_FP_MUL) que mulps. Mais en scalaire, seulement la première sous-unité est sollicitée.


Je vais balancer du code pour illustrer ça, mais là j'ai pas trop le temps

[vectra]

Toujours en apparté, je viens de trouver une lib' pour implémenter des fonctions transcendantales en SSE.
C'est marrant, parce que je me posais justement la question pour accélérer des filtres 2-D complémentaires sur lesquels je travaille là (mais sans le temps de les optimiser).

http://gruntthepeon.free.fr/ssemath/

Il parle aussi d'un type qui bosse sur une version plus rapide encore de log en flottants.
Je verrai ça selon, mais ça m'intéresse dans l'absolu de savoir comment c'est implémenté en standard et en rapide.

[Møgluglu]

Ça a le mérite de tenir dans un seul fichier et d'être KISS, ce qui est pas mal. Après les algos datent des années 80, mais on y revient de toute façon.

Sinon tu as Yeppp!, la bibliothèque dont la mascotte illustre à merveille la légèreté de son marketing :
http://www.yeppp.info/

[vectra]

J'aime pas trop ce type de présentation, on dirait un site orienté développement Android
Mais justement, j'ai l'impression que les contraintes d'embarqué des petits terminaux contribuent nettement à la généralisation des SSE.

La bibliothèque a l'air vraiment intéressante sinon, je vais regarder ça et probablement inclure des papiers du mec. Clairement, on dirait pas que c'est un résultat de recherche...


Sinon, je peux dire dans ma thèse que X86 a largement contribué à en écrire toute une section?

[vectra]

Bon, j'ai toujours pas trouvé le temps de fignoler une présentation de code SIMD pour les nouveaux, mais j'ai déjà des questions pour les grands anciens


Voici une petite fonction toute simple. En entrée, on a un tableau de doubles (r, i) à la suite, et il y a length paires.
Dans le cadre d'un caclul de Fourier, on doit normaliser ces coefficients par un facteur fixe (length, soit le nombre de paires).

J'ai testé en scalaire bête et méchant, en utilisant le type fftw_complex pour représenter le tableau. Je fais 1000 traitements en 0.10s.
Je suis passé en SSE 128 bits, donc un nombre complexe (r, i) par vecteur. Ca déraille total, on passe à 0.40s environ. Et ce, alors que je streame les lectures...
Je suis passé en AVX 256 bits, sur mon PC personnel qui est un Sandy Bridge 2500K (au taf, j'ai du Nehalem).

J'ai donc codé ça:

void
my_local_normalize256(double* __restrict__ data_d, size_t length)
{
  __m256d one = _mm256_set1_pd(double(1.0f));
  __m256d len = _mm256_set1_pd(double(length));
  __m256d div_coef = _mm256_div_pd(one, len);

  size_t run_size = length * 2; // le tableau est complexe
  
#pragma omp parallel for num_threads( 3 )
  for (size_t i = 0; i < run_size; i+=4)
    {
      __m256d cx1 = _mm256_load_pd( &data_d[i] );   
      __m256d res1 = _mm256_mul_pd(cx1, div_coef);
      _mm256_stream_pd(&data_d[i], res1);
    }
#pragma omp barrier
}

Ca plante complet si je prends un tableau de données alloué par fftw, je pense que c'est dû au fait qu'il n'aligne pas sur 32 bits les adresses. J'ai testé avec les lectures-écritures non alignées et non streamées et ça marche. J'ai choisi de remplacer l'allocation de fftw par une simple allocation de tableau _mm_malloc et alignement à 32 bits, et effectivement ça fonctionne.

Mais mal: il me faut 0.40s pour 1000 passes

J'ai donc écrit et testé la même version en non-alignée, non streamée, et là j'ai 0.10s, soit exactement ce que donne le code scalaire.

J'ai donc voulu dumper le code en assembleur pour comprendre ce qui se passe. Petit souci avec la version scalaire: godbolt ne reconnait pas fftw.
Y'a un moyen de "beautifier" du code ASM produit pas gcc? J'aime bien la sortie à godbolt, mais celle de gcc est encore un peu trop roots pour moi.

[Møgluglu]

Si j'ai bien compris, tu as une version avec non-temporal stores qui est 4 fois plus lente que la version avec stores normaux ?
Vu que tu lis et écris le même tableau, ça semble compréhensible.

Je suppose que tout est aligné. Au premier load ou même avant (prefetch), le core va remonter une ligne de cache de 64 octets dans le L1, et te donner les premiers 32 octets dans un registre AVX.

Tu calcules un truc avec, puis tu fais un non-temporal store à la même adresse. Là tu dis au CPU que tu ne vas pas réutiliser ces données, donc il t'obéit servilement en dégageant la ligne concernée des caches pour l'envoyer sur un buffer de write-combining.

À l'itération d'après, tu fais un load sur l'autre moitié de la même ligne de cache, contrairement à ce que tu avais promis. Bim, tu te bouffes un cache miss dans les dents.

tl;dr: les streaming stores c'est seulement pour le streaming, faut pas en mettre pour des algos qui font du read-modify-write en place.

[vectra]

Oui, en effet, tout est aligné.

Eh, c'est pile ça.
C'est à peu près la seule fois que je peux écraser les données d'origine, je n'y avais simplement pas pensé.

Bref, je peux rester sur du load et store alignés classiques, et tenter un déroulement de gros porc (peux pas tester là, je suis sur Nehalem au boulot).


En tous cas, ce fut très pédagogique comme expérience
C'est bien que ça explose quand on fait exploser
Je comprends en cascade toutes les conneries qui sont arrivées, y compris sur le SSE streamé.

Y'a quand-même une raison pour laquelle je voulais streamer, c'est que je vais y être forcé dans l'étape suivante.
En gros, j'utilise mm_shuffle pour séparer les composantes dans deux tableaux distincts R et I, et donc là, stream obligatoire. J'avais dans l'idée de cumuler la séparation streamée et la normalisation dans le même parcours, vu que ça marche dans d'autres cas. Et donc, je me suis acharné à contre-emploi sur ce cas idiot

[vectra]

En revenant sur la différence entre streamé et non-streamé, j'ai constaté qu'un code scalaire de normalisation allait plus vite lorsqu'on lit et écrit au même endroit.

Dans ce cas de figure, on doit de toute façon lire la ligne de cache en entrée pour obtenir les données à normaliser, donc au moment d'écrire cette ligne de cache, on n'aura pas à la relire. Si on était en streaming (entrée et sorties différentes), le non-temporal write éviterait alors cette lecture inutile avant l'écriture.

Mais là, je suppose que les perfs en scalaires in-place doivent être les mêmes que pour un vectoriel non-temporel out-of-place, vu que ce sont les mêmes quantités demouvements de ligne de cache qui sont alors effectueés (par le cache en scalaire, par le cache et le WC en vectoriel)

Sinon, j'ai sorti un code assembleur lisible de la version scalaire / fftw. J'ai l'impression que le code produit est complètement vectorisé et à peu près imbattable:

.LFB12:
	.cfi_startproc
	testq	%rsi, %rsi	# length
	js	.L2	#,
	pxor	%xmm0, %xmm0	# D.4233
	cvtsi2ssq	%rsi, %xmm0	# length, D.4233
.L3:
	testq	%rsi, %rsi	# length
	movss	.LC0(%rip), %xmm1	#, tmp110
	divss	%xmm0, %xmm1	# D.4233, D.4233
	cvtss2sd	%xmm1, %xmm1	# D.4233, dim_factor_inv
	je	.L1	#,
	unpcklpd	%xmm1, %xmm1	# vect_cst_.8
	salq	$4, %rsi	#, D.4234
	xorl	%eax, %eax	# ivtmp.15
	.p2align 4,,10
	.p2align 3
.L5:
	movupd	(%rdi,%rax), %xmm0	# MEM[base: data_8(D), index: ivtmp.15_22, offset: 0B], vect__11.6
	mulpd	%xmm1, %xmm0	# vect_cst_.8, vect__12.7
	movups	%xmm0, (%rdi,%rax)	# vect__12.7, MEM[base: data_8(D), index: ivtmp.15_22, offset: 0B]
	addq	$16, %rax	#, ivtmp.15
	cmpq	%rsi, %rax	# D.4234, ivtmp.15
	jne	.L5	#,
.L1:
	rep ret
	.p2align 4,,10
	.p2align 3
.L2:
	movq	%rsi, %rax	# length, tmp107
	movq	%rsi, %rdx	# length, tmp108
	pxor	%xmm0, %xmm0	# tmp106
	shrq	%rax	# tmp107
	andl	$1, %edx	#, tmp108
	orq	%rdx, %rax	# tmp108, tmp107
	cvtsi2ssq	%rax, %xmm0	# tmp107, tmp106
	addss	%xmm0, %xmm0	# tmp106, D.4233
	jmp	.L3	#
	.cfi_endproc

Faudrait juste remplacer les movu par mova, mais par contre je comprends pas pourquoi l'écriture en mémoire est effectuée par un movupS et pas pD.

[Thamior]

Fondamentalement c'est la même instruction non? Vu que c'est packed le fait d'avoir 2 ou 4 éléments ne te joue pas de tour au niveau du little endian (i.e., recharger derrière avec movupd ou movups c'est pareil si je ne dis pas de bêtises).
Apparemment movups fait 1 byte de moins (internet l'a dit), même si la guideline d'Intel c'est d'utiliser le type qui correspond aux "vraies" données.

De toute façon l'optimisation c'est de la magie noire, tu t'attendais à quoi ?

[vectra]

Je m'attendais clairement à un lock topic si je postais de l'AVX suivi de l'ASM sur le topic de la programmation, déjà
Pourquoi pas une sortie IACA en mode texte, tant qu'on y est?

[Møgluglu]

Apparemment movups fait 1 byte de moins

J'adore.

Quant à mova vs movu, depuis Nehalem ça change que dalle quand les données sont alignées.
Et quand elles ne sont pas alignées movu perd légèrement en perfs tandis que mova crashe. Autant utiliser movu tout le temps (quoique l'instruction fait peut-être un byte du plus ).

[vectra]

Mais comme ça, on peut déplier le code encore un peu plus avant de saturer le cache d'instructions?

[vectra]

Yop les canards noirs. J'ai encore des questions

Il y a un truc que je comprends pas avec LOAD/STORE.

movuss, movups, movaps, ce sont des instructions implémentées par les unités LOAD et STORE qui sont présentes dans le backend CPU



J'ai vu que le nombre de ces unités avait augmenté entre Nehalem et Sandy Bridge. Les commentateurs expliquent que ça améliore les accès-mémoire. J'ai constaté que, en calcul mono-thread, l'impact était clair et net, mais que sur 2 ou 3 threads (en dual-channel), on retombait à peu près dans les mêmes mesures que sur Nehalem.

Je comprends pas bien le rôle de ces unités. Elles sont supposées lire des données dans le cache L1 et les mettre dans des registres? Mais dans ce cas, à quoi ça sert d'en avoir plusieurs par coeur, vu qu'il ne pourra exécuter qu'un seul movaps à la fois?

Au passage, il y en a des tonnes sur les GPU (genre 16 par Streaming Multiprocessor de 32 coeurs).

[Møgluglu]

Je comprends pas bien le rôle de ces unités. Elles sont supposées lire des données dans le cache L1 et les mettre dans des registres? Mais dans ce cas, à quoi ça sert d'en avoir plusieurs par coeur, vu qu'il ne pourra exécuter qu'un seul movaps à la fois?

Tu peux faire plusieurs mov à la fois, si, au moins en scalaire/SSE (il y a une histoire qu'AVX utilise les deux ports de 128 bits à la fois avant Haswell, mais c'est seulement AVX). D'où les 48 bytes/cycle de bande passante du L1. Dans du code x86 typique tu as plein plein de loads depuis la pile d'appel, donc pouvoir traiter 2 loads par cycle depuis le L1 est très utile (et très coûteux en passant ).

Au passage, il y en a des tonnes sur les GPU (genre 16 par Streaming Multiprocessor de 32 coeurs).

Attention, ils comptent chaque voie SIMD comme une unité : c'est plutôt une seule unité de 16 de large. Par contre elle gère les gather/scatter bien plus efficacement qu'un vulgaire Sandy Bridge.

[vectra]

Ah d'acc.

Donc, le L1 est toujours nourri aussi lentement par le prefetch, selon la bande-passante du bus mémoire.
Mais le nombre de load/store va permettre de manger la ligne de cache L1 toute fraichement arrivée beaucoup plus vite. Ca va partir plus vite dans des registres, et donc ça va permettre aux unités de traitement de bosser plus tôt.

Ce qui semble expliquer la super-perf à un coeur, qui s'assagit à 2-3 coeurs.


Un CPU à 3.5 Ghz avec une mémoire à 25 Go/s, le contrôleur mémoire débite au plus 58 bits par cycle de toute manière.
Là, les 48 octets/cycle du L1, c'est sa bande-passante pure pour qui peut le lire à cette vitesse (dans les 150 go/s environ).

Au fait, c'est un cycle, le load/store?

[Møgluglu]

C'est surtout le principe d'un cache. À la base il sert à garder les données qui sont accédées souvent. Si tu as un taux de hit de 90% dans le L1, tu vas utiliser 10x plus de bande passante entre les unités d'exécution et le L1 qu'entre le L1 et le L2.

Ton code qui fait du streaming est une exception, la plupart des applis ont (heureusement) bien plus de localité et profitent des caches. La différence de perf que tu observes sur un seul cœur entre Nehalem et Sandy Bridge doit plutôt venir du débit entre L2 et L3.

Un load/store, c'est 3 ou 4 cycles, mais tu as du store-load forwarding qui te court-circuite le chemin quand tu relis une valeur que tu viens d'écrire, quand les planètes sont bien alignées.

[vectra]

Vache, je manque vraiment de doc là dessus.

Une ligne de cache, c'est toujours 64 octets en L1, L2 et L3?

[Thamior]

Yes sir.

Ah et sinon le store-load forwarding c'est pas gratos non plus .

[Møgluglu]

Ah, le blog d'Henry Wong. Le gars qui arrive à se faire un nom en micro-archi GPU avec deux projets d'étudiant, puis qui il a vu la lumière et qui va faire sa thèse sur du FPGA avec les stars du domaine, tout en continuant à faire de la micro-archi superscalaire et de la programmation d'OS juste pour le fun.

[vectra]

qui il a vu la lumière et qui va faire sa thèse sur du FPGA

Je sens une pointe d'ironie, non?

[Møgluglu]

Non, juste de jalousie, c'est cool les FPGA.

Mais j'ai hésité à ajouter "encore un qui est pas prêt de soutenir".