Accélérez la multiplication de la matrice flottante 4x4 avec SIMD

De nombreuses années se sont écoulées depuis que je me suis familiarisé avec les instructions MMX, SSE et plus tard AVX sur les processeurs Intel. À un moment donné, ils semblaient être une sorte de magie dans le contexte de l'assembleur x86, qui avait longtemps été quelque chose de banal. Ils m'ont tellement accroché qu'il y a quelques années, j'ai eu l'idée d'écrire mon propre logiciel de rendu pour un jeu célèbre. La chose qui a promis ces instructions m'a promis ceci. À un moment donné, j'ai même pensé à l'écrire. Mais l'écriture de texte s'est avérée beaucoup plus compliquée que le code.

À cette époque, je voulais éviter les problèmes de support sur différents processeurs. Je voulais pouvoir vérifier mon rendu sur le montant maximum disponible. J'ai toujours des amis avec d'anciens processeurs AMD, et leur plafond était SSE3. Par conséquent, à ce moment-là, j'ai décidé de me limiter à un maximum de SSE3. Il y avait donc une bibliothèque mathématique vectorielle, un peu moins complètement implémentée sur SSE, avec une inclusion rare avant SSE3. Cependant, à un moment donné, je me suis demandé quelles performances maximales je pouvais retirer du processeur pour un certain nombre d'opérations mathématiques vectorielles critiques. Une telle opération consiste à multiplier 4 matrices flottantes par 4.


En fait, j'ai décidé de faire cette affaire davantage pour le plaisir. J'ai déjà écrit et j'utilise la multiplication matricielle pour mon rendu logiciel sur SSE et cela me semble suffisant. Mais j'ai alors décidé de voir combien de mesures je pouvais retirer en principe en multipliant 2 matrices float4x4. Sur mon SSE actuel, ce sont 16 cycles d'horloge. Certes, la transition récente vers IACA 3 a commencé à afficher 19, car j'ai commencé à écrire 1 * pour certaines instructions au lieu de 0 *. Apparemment plus tôt, ce n'était qu'un défaut dans l'analyseur.

En bref sur les utilitaires utilisés

Pour l'analyse de code, j'ai utilisé le célèbre utilitaire Intel Architecture Code Analyzer . Pour l'analyse, j'utilise l'architecture Haswell (HSW), au minimum avec le support d'AVX2. Pour écrire du code est également très pratique à utiliser: Intel Intrinsics Guide et Intel optimization manual .

Pour l'assemblage, j'utilise la communauté MSVS 2017 à partir de la console. J'écris le code dans la version avec intrinsèques. Vous écrivez une fois, et généralement cela fonctionne immédiatement sur différentes plateformes. De plus, le compilateur x64 VC ++ ne prend pas en charge l'assembleur en ligne, mais je veux qu'il fonctionne également sous x64.

Étant donné que cet article est déjà un peu au-delà du niveau débutant en programmation SIMD, je ne décrirai pas les registres, les instructions, dessiner (ou raser) de belles images et essayer d'apprendre la programmation en utilisant les instructions SIMD. Le site Web d'Intel regorge d'une documentation excellente, claire et détaillée.

Je voulais tout rendre plus facile ... Mais il s'est avéré comme toujours

C'est là que le moment commence, ce qui complique beaucoup la mise en œuvre et l'article. Je m'y attarderai donc un peu. Ce n'est pas intéressant pour moi d'écrire une multiplication matricielle avec une disposition standard des éléments en ligne. Qui en avait besoin, et donc ils ont étudié dans les universités ou seuls. Notre objectif est la productivité. Tout d'abord, je suis passé à la disposition des colonnes il y a longtemps. Mon logiciel de rendu est basé sur l'API OpenGL et donc, afin d'éviter des transpositions inutiles, j'ai commencé à stocker des éléments dans des colonnes. Ceci est également important car la multiplication matricielle n'est pas si critique. Matrices 2-5-10 bien multipliées. Et c'est tout. Et puis nous multiplions la matrice finie par des milliers ou des millions de sommets. Et cette opération est beaucoup plus critique. Vous pouvez bien sûr transposer à chaque fois. Mais pourquoi, si cela peut être évité.

Mais revenons aux matrices exclusivement. Nous avons déterminé le stockage en colonnes. Cependant, cela peut être encore plus compliqué. Il est plus pratique pour moi de stocker les éléments supérieurs des vecteurs et des lignes de matrice dans les registres SIMD afin que x soit dans le flottant le plus élevé (index 3) et w dans le mineur (index 0). Ici, apparemment, nous devrons retraiter à nouveau sur pourquoi.

Le fait est que dans un logiciel de rendu dans un vecteur, vous devez manipuler le composant w plus souvent ( 1 / z y est stocké), et il est très pratique de le faire via la version _ss de l'opération (opérations exclusivement avec le composant dans le flotteur inférieur du registre xmm ), sans toucher x, y, z . Par conséquent, dans le registre SSE, le vecteur est stocké dans un ordre compréhensible x, y, z, w et en mémoire dans le sens inverse w, z, y, x .

De plus, toutes les options de multiplication sont également implémentées par des fonctions individuelles. Cela se fait parce que je les utilise pour remplacer l'option souhaitée en fonction du type d'instructions prises en charge. Bien décrit ici.

Nous implémentons la fonctionnalité de base

Multiplication avec boucles, rangée ordonnée

Ici, tout est simple et clair. Pour chaque élément, nous faisons 4 multiplications et 3 additions. Au total, ce sont 64 multiplications et 48 additions. Et cela sans tenir compte de la lecture des éléments d'enregistrement.

Tout est triste, bref. Pour cette option, pour le cycle interne, l'IACA a émis: 3,65 cycles d'horloge pour l'assemblage x86 et 2,97 horloges pour l'assemblage x64 . Ne demandez pas pourquoi les nombres fractionnaires. Je ne sais pas. IACA 2.1 n'en a pas souffert. Dans tous les cas, ces nombres doivent être multipliés par environ 4 * 4 * 4 = 64. Même si vous prenez x64, le résultat est d'environ 192 mesures. Il est clair qu'il s'agit d'une estimation approximative. Je ne vois pas l'intérêt d'évaluer plus précisément les performances de cette option.

Implémentation en boucle, colonne ordonnée

Multiplication de cycle, stockage orienté SIMD

Cette implémentation simplifie quelque peu la compréhension de ce qui se passe à l'intérieur, mais n'est clairement pas suffisante.

Classes d'assistance

Pour faciliter la compréhension et l'écriture de code de référence et de débogage, il est pratique d'implémenter quelques classes auxiliaires. Rien de plus, tout n'est que pour comprendre. Je note que la mise en œuvre de classes vectorielles et matricielles à part entière est une question difficile distincte, et n'est pas incluse dans le sujet de cet article.

Fonction de référence pour les tests

Étant donné que l'ordre accepté des éléments dans la matrice complique beaucoup la compréhension, nous ne serons pas non plus dérangés par la fonction de référence claire , qui montrera dans les futures implémentations que tout fonctionne correctement. Nous comparerons avec lui les résultats ultérieurs.


L'algorithme classique est clairement peint ici, il est difficile de se tromper (mais vous pouvez :-)). À ce sujet, l'IACA a publié: x86 - 69,95 mesures, x64 - 64 mesures . Voici environ 64 cycles et nous verrons l'accélération de cette opération dans le futur.

Implémentation SSE

Algorithme SSE classique

Pourquoi classique? Parce qu'il fait depuis longtemps partie de l'implémentation de FVec dans le cadre de MSVS. Pour commencer, nous allons écrire comment nous présentons les éléments de matrice dans les registres SSE. Cela semble déjà plus simple ici. Juste une matrice transposée.

 //     00, 10, 20, 30 // m[0] -  SIMD /   01, 11, 21, 31 // m[1] 02, 12, 22, 32 // m[2] 03, 13, 23, 33 // m[3] 

Nous prenons le code de déroulement de la variante ci-dessus. D'une certaine manière, il est hostile à SSE. Le premier groupe de lignes comprend les résultats de la colonne de la matrice résultante: r._00, r._01, r._02, r._03 . Nous avons cette colonne, mais nous avons besoin d'une ligne. Oui, et m , n ne semblent pas pratiques pour les calculs. Par conséquent, nous réorganisons les lignes de l'algorithme de sorte que le résultat r soit rangé.

 //  ,     r[0] r00 = m00*n00 + m01*n10 + m02*n20 + m03*n30; r10 = m10*n00 + m11*n10 + m12*n20 + m13*n30; r20 = m20*n00 + m21*n10 + m22*n20 + m23*n30; r30 = m30*n00 + m31*n10 + m32*n20 + m33*n30; //  ,     r[1] r01 = m00*n01 + m01*n11 + m02*n21 + m03*n31; r11 = m10*n01 + m11*n11 + m12*n21 + m13*n31; r21 = m20*n01 + m21*n11 + m22*n21 + m23*n31; r31 = m30*n01 + m31*n11 + m32*n21 + m33*n31; //  ,     r[2] r02 = m00*n02 + m01*n12 + m02*n22 + m03*n32; r12 = m10*n02 + m11*n12 + m12*n22 + m13*n32; r22 = m20*n02 + m21*n12 + m22*n22 + m23*n32; r32 = m30*n02 + m31*n12 + m32*n22 + m33*n32; //  ,     r[3] r03 = m00*n03 + m01*n13 + m02*n23 + m03*n33; r13 = m10*n03 + m11*n13 + m12*n23 + m13*n33; r23 = m20*n03 + m21*n13 + m22*n23 + m23*n33; r33 = m30*n03 + m31*n13 + m32*n23 + m33*n33; 

Mais c'est déjà beaucoup mieux. Que voit-on en fait? Selon les colonnes de l'algorithme dans chaque groupe, nous avons les lignes de la matrice m impliquées:

 m [0] = {00,10,20,30}, m [1] = {01,11,21,31}, m [2] = {02,12,22,32}, m [3] = {03,13,23,33},

qui sont multipliés par le même élément de la matrice n . Par exemple, pour le premier groupe, c'est: n._00, n._10, n._20, n._30 . Et les éléments de la matrice n pour chaque groupe de lignes de l'algorithme se trouvent à nouveau dans une ligne de la matrice.

Alors tout est simple: on prend simplement les lignes de la matrice m par index, mais comme pour les éléments n , on prend sa ligne et on la mélange à travers les 4 éléments du registre via l'instruction shuffle , pour multiplier par la ligne de la matrice m dans le registre. Par exemple, pour l'élément n._00 (rappelez-vous que son décalage dans le registre a l'index 3), ce sera:

  _mm_shuffle_ps (n [0], n [0], _MM_SHUFFLE (3,3,3,3)) 

Dans une forme simplifiée, l'algorithme ressemble à ceci:

 //   n[0]={00,10,20,30} r[0] = m[0] * n00 + m[1] * n10 + m[2] * n20 + m[3] * n30; //   n[1]={01,11,21,31} r[1] = m[0] * n01 + m[1] * n11 + m[2] * n21 + m[3] * n31; //   n[2]={02,12,22,32} r[2] = m[0] * n02 + m[1] * n12 + m[2] * n22 + m[3] * n32; //   n[3]={03,13,23,33} r[3] = m[0] * n03 + m[1] * n13 + m[2] * n23 + m[3] * n33; 

Maintenant, nous changeons les éléments n de l'algorithme en shuffle correspondant, multiplication par _mm_mul_ps , la somme par _mm_add_ps , et vous avez terminé. Ça marche. Le code semble cependant bien pire que l'algorithme lui-même. Pour ce code, l'IACA a émis: x86 - 18.89, x64 - 16 cycles . C'est 4 fois plus rapide que le précédent. Dans le registre SSE, le quatrième flotteur. Relation presque linéaire.

Décorer la mise en œuvre de SSE

Pourtant, dans le code, cela a l'air horrible. Nous allons essayer d'améliorer cela en écrivant un peu de sucre syntaxique.


Le compilateur peut parfaitement intégrer ces fonctions (bien que parfois sans __forceinline en aucune façon).


Et donc c'est déjà beaucoup mieux et plus lisible. Pour cela, IACA a produit environ le résultat attendu: x86 - 19 (et pourquoi pas fractionnaire?), X64 - 16 . En fait, les performances n'ont pas changé, mais le code est beaucoup plus beau et compréhensible.

Peu de contribution à l'optimisation future

Introduisons encore une amélioration au niveau d'une fonction apparue récemment dans la version iron. Opération d'ajout multiple (fma) . fma (a, b, c) = a * b + c .


Pourquoi est-ce nécessaire? Tout d'abord, pour une optimisation future. Par exemple, vous pouvez simplement remplacer mad par fma dans le code fini via les mêmes macros que vous le souhaitez. Mais nous allons jeter les bases de l'optimisation maintenant:


IACA: x86 - 18,89, x64 - 16 . Encore une fois fractionnaire. Pourtant, IACA produit parfois des résultats étranges. Le code n'a pas tellement changé. Probablement même un peu pire. Mais l'optimisation nécessite parfois de tels sacrifices.

Nous passons à l'épargne via _mm_stream

Divers guides d'optimisation recommandent à nouveau de ne pas extraire le cache pour les opérations de sauvegarde en bloc. Cela est généralement justifié lorsque vous traitez des sommets qui sont des milliers ou plus. Mais pour les matrices, ce n'est peut-être pas si important. Cependant, je vais l'ajouter quand même.


Rien n'a changé dans le temps ici, du mot du tout. Mais, selon les recommandations, nous ne touchons plus au cache.

Implémentation AVX

Option AVX de base

Ensuite, nous passons à la prochaine étape d'optimisation. Le 4ème flotteur est inclus dans le registre SSE, et dans l'AVX il est déjà 8. Autrement dit, il y a une chance théorique de réduire le nombre d'opérations effectuées et d'augmenter la productivité sinon de moitié, puis au moins 1,5 fois. Mais quelque chose me dit que tout ne sera pas aussi simple avec la transition vers AVX. Pouvons-nous obtenir les données nécessaires à partir de registres doubles?

Essayons de le comprendre. Encore une fois, nous écrivons notre algorithme de multiplication utilisé ci-dessus. Vous ne pouvez pas faire cela, mais il est plus pratique de traiter le code lorsque tout est à proximité et que vous n'avez pas à faire défiler une demi-page.

 //    : 00, 10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33 //   SSE: r0 = m0*n00 + m1*n10 + m2*n20 + m3*n30 r1 = m0*n01 + m1*n11 + m2*n21 + m3*n31 r2 = m0*n02 + m1*n12 + m2*n22 + m3*n32 r3 = m0*n03 + m1*n13 + m2*n23 + m3*n33 

En sortie, nous nous attendons à obtenir le résultat en ymm = {r0: r1} et ymm = {r2: r3} . Si dans la version SSE notre algorithme était généralisé aux colonnes, nous devons maintenant le généraliser aux lignes. Donc, agir comme dans le cas de l'option SSE ne fonctionnera pas.

Si l'on considère la matrice m dans les registres ymm , on obtient respectivement ymm = {m0: m1} et ymm = {m2: m3} . Auparavant, nous n'avions que des colonnes matricielles dans le registre, et maintenant des colonnes et des lignes.

Si vous essayez d'agir comme auparavant, vous devez multiplier ymm = {m0: m1} par le registre ymm = {n00, n00, n00, n00}: {n10, n10, n10, n10} . Étant donné que n00 et n01 sont dans la même ligne de la matrice n , à en juger par l'ensemble d'instructions AVX disponibles, les disperser par ymm sera coûteux. La lecture aléatoire et la permutation fonctionnent séparément pour chacun des deux fours d'un flotteur (haut et bas xmm ) à l'intérieur des registres ymm .

Si nous prenons ymm de la matrice n , alors nous obtenons les deux éléments n00 et n10 au plus haut de 2 xmm à l'intérieur du registre ymm . {n00, n10, n20, n30}: {n01, n11, n21, n31} . Habituellement, l'index pour les instructions existantes est de 0 à 3. Et il ne flotte que dans un registre xmm sur deux dans le registre ymm . Il n'est pas possible de transférer le n10 de l'ancien xmm vers le plus jeune à moindre coût . Et puis cet objectif doit être répété plusieurs fois. Nous ne pouvons pas supporter une telle perte de mesures. Il faut trouver autre chose.

Nous généralisions les colonnes, mais maintenant les lignes. Par conséquent, nous essaierons d'aller d'une manière un peu différente. Nous devons obtenir le résultat dans {r0: r1} . Cela signifie que l'algorithme doit être amélioré non pas sur des lignes distinctes de l'algorithme, mais sur deux à la fois. Et ici, ce qui était un inconvénient dans le travail de mélange et de permutation , sera un plus pour nous. Nous regardons ce que nous aurons dans les registres ymm lorsque nous considérons la matrice n .

 n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31} n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33} 

Oui, nous remarquons que dans différentes parties xmm du registre ymm, nous avons les éléments 00 et 01 . Ils peuvent être multipliés en casse par registre via la commande permute dans {_00, _00, _00, _00}: {_ 01, _01, _01, _01} , indiquant un seul index 3 pour les deux parties xmm . C'est exactement ce dont nous avons besoin. En effet, les coefficients sont également utilisés sur différentes lignes. Seulement maintenant dans le registre ymm correspondant pour la multiplication, il sera nécessaire de conserver {m0: m0} , c'est-à-dire la première ligne dupliquée de la matrice m .

Donc, nous peignons l'algorithme plus en détail. Nous lisons les doubles rangées de la matrice m dans les registres ymm :

 mm[0] = {m0:m0} mm[1] = {m1:m1} mm[2] = {m2:m2} mm[3] = {m3:m3} 

Et puis nous calculerons la multiplication comme:

 r0r1 = mm[0] * {n00,n00,n00,n00:n01,n01,n01,n01} + // permute<3,3,3,3>(n0n1) mm[1] * {n10,n10,n10,n10:n11,n11,n11,n11} + // permute<2,2,2,2>(n0n1) mm[2] * {n20,n20,n20,n20:n21,n21,n21,n21} + // permute<1,1,1,1>(n0n1) mm[3] * {n30,n30,n30,n30:n31,n31,n31,n31} // permute<0,0,0,0>(n0n1) r2r3 = mm[0] * {n02,n02,n02,n02:n03,n03,n03,n03} + // permute<3,3,3,3>(n2n3) mm[1] * {n12,n12,n12,n12:n13,n13,n13,n13} + // permute<2,2,2,2>(n2n3) mm[2] * {n22,n22,n22,n22:n23,n23,n23,n23} + // permute<1,1,1,1>(n2n3) mm[3] * {n32,n32,n32,n32:n33,n33,n33,n33} // permute<0,0,0,0>(n2n3) 

Nous réécrivons plus clairement:

 r0r1 = mm[0]*n0n1<3,3,3,3>+mm[1]*n0n1<2,2,2,2>+mm[2]*n0n1<1,1,1,1>+mm[3]*n0n1<0,0,0,0> r2r3 = mm[0]*n2n3<3,3,3,3>+mm[1]*n2n3<2,2,2,2>+mm[2]*n2n3<1,1,1,1>+mm[3]*n2n3<0,0,0,0> 

Ou sous une forme simplifiée:

 r0r1 = mm[0]*n0n1<3> + mm[1]*n0n1<2> + mm[2]*n0n1<1> + mm[3]*n0n1<0> r2r3 = mm[0]*n2n3<3> + mm[1]*n2n3<2> + mm[2]*n2n3<1> + mm[3]*n2n3<0> 

Tout semble clair.


Voici les chiffres intéressants de l'IACA: x86 - 12,53, x64 - 12 . Bien sûr, je voulais mieux. J'ai raté quelque chose.

Optimisation AVX plus sucre syntaxique

Il semble que dans le code ci-dessus, AVX n'ait pas été utilisé à son plein potentiel. Nous constatons qu'au lieu de définir deux lignes identiques dans le registre ymm , nous pouvons utiliser la diffusion , qui peut remplir le registre ymm avec deux valeurs xmm identiques. En cours de route, ajoutez également du «sucre syntaxique» pour les fonctions AVX.


Et ici, les résultats sont déjà plus intéressants. IACA produit des nombres: x86 - 10, x64 - 8,58 , ce qui semble beaucoup mieux, mais toujours pas 2 fois.

Option AVX + FMA (finale)

Faisons une autre tentative. Il serait maintenant logique de rappeler à nouveau le jeu d'instructions FMA, car il a été ajouté aux processeurs après AVX. Il suffit de changer individuellement mul + add pour une seule opération. Bien que nous utilisions toujours l'instruction de multiplication pour donner au compilateur plus de possibilités d'optimisation et le processeur pour l'exécution parallèle de multiplications. Habituellement, je regarde le code généré dans l'assembleur pour m'assurer quelle option est meilleure.

Dans ce cas, nous devons calculer a * b + c * d + e * f + g * h . Vous pouvez faire ce front: fma (a, b, fma (c, d, fma (e, f, g * h))) . Mais, comme nous le voyons, il est impossible d'effectuer une opération ici sans terminer la précédente. Et cela signifie que nous ne pourrons pas utiliser la possibilité de faire des multiplications par paires, comme le permet le pipeline SIMD. Si nous transformons légèrement les calculs fma (a, b, c * d) + fma (e, f, g * h) , nous verrons que nous pouvons paralléliser les calculs. Faites d'abord deux multiplications indépendantes, puis deux opérations fma indépendantes.


IACA: x86 - 9.21, x64 - 8 . Maintenant c'est très bien. Quelqu'un dira probablement ce qui peut être fait encore mieux, mais je ne sais pas comment.

Repères

Je constate tout de suite que ces chiffres ne doivent pas être considérés comme la vérité ultime. Même avec un test fixe, ils nagent dans certaines limites.Et plus encore, ils se comportent différemment sur différentes plateformes. Avec toute optimisation, prenez des mesures spécifiquement pour votre cas.

Table des matières

• Fonction: nom de la fonction. La fin sur s - fonctionne avec un streaming, mov autrement normal (sans streaming). Ajouté pour plus de clarté, car c'est assez important.
• Cycles IACA: nombre de ticks par fonction calculé par IACA
• Cycles mesurés: nombre mesuré de mesures (moins c'est plus)
• Accélération IACA: nombre de mesures sur une ligne zéro / nombre de mesures sur une ligne
• Accélération mesurée: nombre de mesures dans la ligne zéro / nombre de mesures dans la ligne (plus c'est mieux)

Pour loop_m, les graduations de l'article ont été multipliées par 64. Autrement dit, il s'agit d'une valeur très approximative. En fait, cela s'est avéré ainsi.

i3-3770:

i7-8700K:

Testez le code dans la source. S'il existe des suggestions raisonnables sur la façon de les améliorer, écrivez dans les commentaires.

BONUS du royaume de la fiction

En fait, c'est du domaine de la fiction parce que si je voyais des processeurs prenant en charge AVX512, alors peut-être sur les photos. Cependant, j'ai essayé d'implémenter l'algorithme. Ici, je ne vais rien expliquer, une analogie complète avec AVX + FMA. L'algorithme est le même, seulement moins d'opérations.


Les chiffres sont fantastiques: x86 - 4,79, x64 - 5,42 (IACA avec architecture SKX). Ceci malgré le fait que l'algorithme a 64 multiplications et 48 additions.

Code PS de l'article

Ceci est ma première expérience en écrivant un article. Merci à tous pour vos commentaires. Ils aident à améliorer le code et l'article.