Multiplication masortingcielle 2×2 optimisée: Assemblage lent versus rapide SIMD

Problème

J’étudie des algorithmes de multiplication de masortingce hautes performances tels que OpenBLAS ou GotoBLAS et j’essaie de reproduire certains des résultats. Cette question concerne le kernel interne d’un algorithme de multiplication de masortingces. Plus précisément, je cherche à calculer C += AB , où A et B sont des masortingces 2×2 de type double à la vitesse maximale de mon processeur. Il y a deux façons de faire ça. Une façon consiste à utiliser les instructions SIMD. La deuxième méthode consiste à coder directement dans l’assemblage à l’aide des registres SIMD.

Ce que j’ai regardé jusqu’à présent

Tous les articles pertinents, les pages Web du cours, de nombreuses questions et réponses sur le sujet (trop nombreuses pour être énumérées), j’ai compilé OpenBLAS sur mon ordinateur, parcouru les codes source OpenBLAS, GotoBLAS et BLIS, les manuels d’Agner.

• que fait cette instruction?: – mov% gs: 0x14,% eax
• Comment multiplier deux quaternions avec des instructions minimales?
• AMD64 – Instruction d’assemblage nopw?
• Correction à chaud d’une fonction
• Code d’assemblage fsqrt et instructions fmul

Matériel

Mon processeur est un Intel i5 – 540M. Vous pouvez trouver les informations CPUID pertinentes sur cpu-world.com. La microarchitecture est Nehalem (westmere), elle permet donc théoriquement de calculer 4 bascules à double précision par cœur et par cycle. Je vais utiliser un seul cœur (pas d’OpenMP), donc avec l’hyperthreading et Intel Turbo Boost à 4 étapes, je devrais voir un pic de ( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops . Pour référence, avec les deux cœurs tournant au maximum, Intel Turbo Boost accélère en 2 étapes et je devrais obtenir un maximum théorique de 22.4 DP Gflops .

Installer

Je déclare mes masortingces 2×2 comme double et les initialise avec des entrées aléatoires, comme indiqué dans l’extrait de code ci-dessous.

 srand(time(NULL)); const int n = 2; double A[n*n]; double B[n*n]; double C[n*n]; double T[n*n]; for(int i = 0; i < n*n; i++){ A[i] = (double) rand()/RAND_MAX; B[i] = (double) rand()/RAND_MAX; C[i] = 0.0; } 

Je calcule une vraie réponse en utilisant une multiplication naïve masortingce-masortingce (illustrée ci-dessous) qui me permet de vérifier mon résultat soit visuellement, soit en calculant la norme L2 de tous les éléments.

 // "true" answer for(int i = 0; i < n; i++) for(int j = 0; j < n; j++) for(int k = 0; k < n; k++) T[i*n + j] += A[i*n + k]*B[k*n + j]; 

Pour exécuter le code et obtenir une estimation des Gflops, j’appelle chaque fonction de multiplication une fois pour le réchauffer, puis je l’exécute dans une boucle for pendant un nombre maxiter fois, en veillant à maxiter à zéro la masortingce C chaque fois que je calcule C += AB . La boucle for est placée dans deux instructions clock() et est utilisée pour estimer les Gflops. Le coup de code illustre cette partie.

 C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; mult2by2(A,B,C); //warmup time1 = clock(); for(int i = 0; i < maxiter; i++){ mult2by2(A,B,C); C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; } time2 = clock() - time1; time3 = (double)(time2)/CLOCKS_PER_SEC; gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; mult2by2(A,B,C); // to compute the norm against T norm = L2norm(n,C,T); 

Code SIMD

Mon processeur prend en charge les vecteurs 128 bits, ce qui me permet d’adapter 2 double s dans chaque vecteur. C’est la raison principale pour laquelle je fais la multiplication de masortingce 2×2 dans le kernel interne. Le code SIMD calcule une ligne entière de C à la fois.

  inline void __atsortingbute__ ((gnu_inline)) __atsortingbute__ ((aligned(16))) mult2by2B( const double* ressortingct A, const double* ressortingct B, double* ressortingct C ) { register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; xmm0 = _mm_load_pd(C); xmm1 = _mm_load1_pd(A); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 1); xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C,xmm2); xmm0 = _mm_load_pd(C + 2); xmm1 = _mm_load1_pd(A + 2); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 3); //xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C + 2,xmm2); } 

Assmebly (Intel Syntax)

Ma première tentative a été de créer une routine d’assemblage distincte pour cette pièce et de l’appeler à partir de la routine main . Cependant, cela a été extrêmement lent, car je ne parviens pas à extern fonctions extern . J’ai écrit l’assemblage en tant qu’assemblage en ligne, comme indiqué ci-dessous. Il est identique à celui produit par gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel . D’après ce que j’ai compris des diagrammes de microarchitecture de Nehalem, ce processeur peut exécuter en parallèle SSE ADD , SSE MUL et SSE MOV , ce qui explique l’entrelacement d’instructions MUL , ADD et MOV . Vous remarquerez que les instructions SIMD ci-dessus sont dans un ordre différent, car j’avais une compréhension différente des manuels d’Agner Fog. Néanmoins, gcc est intelligent et le code SIMD ci-dessus est compilé dans l’assemblage indiqué dans la version intégrée.

 inline void __atsortingbute__ ((gnu_inline)) __atsortingbute__ ((aligned(16))) mult2by2A ( const double* ressortingct A, const double* ressortingct B, double* ressortingct C ) { __asm__ __volatile__ ( "mov edx, %[A] \n\t" "mov ecx, %[B] \n\t" "mov eax, %[C] \n\t" "movapd xmm3, XMMWORD PTR [ecx] \n\t" "movapd xmm2, XMMWORD PTR [ecx+16] \n\t" "movddup xmm1, QWORD PTR [edx] \n\t" "mulpd xmm1, xmm3 \n\t" "addpd xmm1, XMMWORD PTR [eax] \n\t" "movddup xmm0, QWORD PTR [edx+8] \n\t" "mulpd xmm0, xmm2 \n\t" "addpd xmm0, xmm1 \n\t" "movapd XMMWORD PTR [eax], xmm0 \n\t" "movddup xmm4, QWORD PTR [edx+16] \n\t" "mulpd xmm4, xmm3 \n\t" "addpd xmm4, XMMWORD PTR [eax+16] \n\t" "movddup xmm5, QWORD PTR [edx+24] \n\t" "mulpd xmm5, xmm2 \n\t" "addpd xmm5, xmm4 \n\t" "movapd XMMWORD PTR [eax+16], xmm5 \n\t" : // no outputs : // inputs [A] "m" (A), [B] "m" (B), [C] "m" (C) : //register clobber "memory", "edx","ecx","eax", "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5" ); } 

Résultats

Je comstack mon code avec les drapeaux suivants:

 gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 

Les résultats pour maxiter = 1000000000 sont ci-dessous:

 ********** Inline ASM L2 norm: 0.000000e+000, Avg. CPU time: 9.563000, Avg. Gflops: 1.673115 ********** SIMD Version L2 norm: 0.000000e+000, Avg. CPU time: 0.359000, Avg. Gflops: 44.568245 

Si je force la version SIMD à ne pas être __atsortingbute__ ((noinline)) avec __atsortingbute__ ((noinline)) , les résultats sont les suivants:

 ********** Inline ASM L2 norm: 0.000000e+000, Avg. CPU time: 11.155000, Avg. Gflops: 1.434334 ********** SIMD Version L2 norm: 0.000000e+000, Avg. CPU time: 11.264000, Avg. Gflops: 1.420455 

Des questions

1. Si les implémentations ASM en ligne et SIMD produisent toutes deux une sortie d’assemblage identique, pourquoi la version d’assemblage est-elle tellement plus lente? C’est comme si l’assemblage en ligne ne s’était pas aligné, ce qui est mis en évidence par le deuxième ensemble de résultats montrant des performances identiques pour l’ASM “en ligne” par rapport au SIMD “sans alignement”. La seule explication que je puisse trouver est dans Agner Fog Volume 2 page 6 :

   Le code compilé peut être plus rapide que le code assembleur car les compilateurs peuvent optimiser les procédures et les programmes dans leur intégralité. Le programmeur d’assemblage doit généralement créer des fonctions bien définies avec une interface d’appel bien définie qui obéit à toutes les conventions d’appel afin de rendre le code testable et vérifiable. Cela empêche bon nombre des méthodes d’optimisation utilisées par les compilateurs, telles que l’inligne de fonction, l’allocation de registre, la propagation constante, l’élimination des sous-expressions communes entre fonctions, la planification entre fonctions, etc. Ces avantages peuvent être obtenus en utilisant du code C ++ avec des fonctions insortingnsèques au lieu du code d’assemblage .

   Mais la sortie de l’assembleur pour les deux versions est exactement la même.

2. Pourquoi vois-je 44 Gflops dans la première série de résultats? C’est bien au-dessus du pic de 12 Gflops que j’ai calculé, et c’est ce que j’attendrais si j’exécutais les deux cœurs avec des calculs à simple précision.

EDIT 1 Le commentaire indique qu’il peut y avoir une élimination du code mort. Je peux confirmer que cela se produit pour les instructions SIMd. La sortie -S montre que la boucle for pour SIMD ne contient que la masortingce C zéros. Je peux le désactiver en désactivant l’optimisation du compilateur avec -O0 . Dans ce cas, SIMD fonctionne 3 fois plus lentement que ASM, mais ASM fonctionne toujours à la même vitesse. La norme est également non nulle à présent, mais ça rest OK à 10 ^ -16. Je vois aussi que la version ASM en ligne est en train d’être alignée avec les balises APP et NO_APP , mais elle est également déroulée 8 fois dans la boucle for . Je pense que le fait de dérouler cela autant de fois aura un impact important sur les performances, comme je le fais habituellement 4 fois. D’après mon expérience, rien de plus semble dégrader les performances.

• Comment utiliser GDB en mode 16 bits?
• Division de bits en cours d’exécution sans opérateurs arithmétiques
• Quelles sont les principales étapes de la compilation?
• Comment l’assemblage affecte-t-il un nombre négatif à un entier non signé?
• Comment le code de démarrage Linux peut-il être écrit en C?

 C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 

 #include  #include  // SSE2 #include  inline void __atsortingbute__ ((gnu_inline)) __atsortingbute__ ((aligned(16))) mult2by2B( const double* __ressortingct A, const double* __ressortingct B, double* __ressortingct C ) { register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; xmm0 = _mm_load_pd(C); xmm1 = _mm_load1_pd(A); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 1); xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C,xmm2); xmm0 = _mm_load_pd(C + 2); xmm1 = _mm_load1_pd(A + 2); xmm2 = _mm_load_pd(B); xmm3 = _mm_load1_pd(A + 3); //xmm4 = _mm_load_pd(B + 2); xmm1 = _mm_mul_pd(xmm1,xmm2); xmm2 = _mm_add_pd(xmm1,xmm0); xmm1 = _mm_mul_pd(xmm3,xmm4); xmm2 = _mm_add_pd(xmm1,xmm2); _mm_store_pd(C + 2,xmm2); } int main() { double A[4], B[4], C[4]; int maxiter = 10000000; //int maxiter = 1000000000; double dtime; dtime = omp_get_wtime(); for(int i = 0; i < maxiter; i++){ mult2by2B(A,B,C); C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; } dtime = omp_get_wtime() - dtime; printf("%f %f %f %f\n", C[0], C[1], C[2], C[3]); //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; printf("time %f\n", dtime); }