Optimisation à l’aide de SSE Insortingnsics

J’essaie de convertir une boucle que j’ai en un élément insortingnsèque du SSE. Je semble avoir fait d’assez bons progrès, et par là, je veux dire que c’est dans la bonne direction, mais il semble que j’ai mal interprété la traduction quelque part, car je n’obtiens pas la même réponse “correcte” qui découle du code non-sse. .

Ma boucle originale que j’ai déroulée par 4 est la suivante:

int unroll_n = (N/4)*4; for (int j = 0; j < unroll_n; j++) { for (int i = 0; i < unroll_n; i+=4) { float rx = x[j] - x[i]; float ry = y[j] - y[i]; float rz = z[j] - z[i]; float r2 = rx*rx + ry*ry + rz*rz + eps; float r2inv = 1.0f / sqrt(r2); float r6inv = r2inv * r2inv * r2inv; float s = m[j] * r6inv; ax[i] += s * rx; ay[i] += s * ry; az[i] += s * rz; //u rx = x[j] - x[i+1]; ry = y[j] - y[i+1]; rz = z[j] - z[i+1]; r2 = rx*rx + ry*ry + rz*rz + eps; r2inv = 1.0f / sqrt(r2); r6inv = r2inv * r2inv * r2inv; s = m[j] * r6inv; ax[i+1] += s * rx; ay[i+1] += s * ry; az[i+1] += s * rz; //unroll i 3 rx = x[j] - x[i+2]; ry = y[j] - y[i+2]; rz = z[j] - z[i+2]; r2 = rx*rx + ry*ry + rz*rz + eps; r2inv = 1.0f / sqrt(r2); r6inv = r2inv * r2inv * r2inv; s = m[j] * r6inv; ax[i+2] += s * rx; ay[i+2] += s * ry; az[i+2] += s * rz; //unroll i 4 rx = x[j] - x[i+3]; ry = y[j] - y[i+3]; rz = z[j] - z[i+3]; r2 = rx*rx + ry*ry + rz*rz + eps; r2inv = 1.0f / sqrt(r2); r6inv = r2inv * r2inv * r2inv; s = m[j] * r6inv; ax[i+3] += s * rx; ay[i+3] += s * ry; az[i+3] += s * rz; } } 

Ensuite, je suis allé ensuite ligne par ligne pour la section supérieure et l’ai convertie en éléments insortingnsèques SSE. Le code est ci-dessous. Je ne suis pas tout à fait sûr que les trois premières lignes soient nécessaires, mais je comprends que mes données doivent être alignées sur 16 bits pour que cela fonctionne correctement et de manière optimale.

 float *x = malloc(sizeof(float) * N); float *y = malloc(sizeof(float) * N); float *z = malloc(sizeof(float) * N); for (int j = 0; j < N; j++) { for (int i = 0; i < N; i+=4) { __m128 xj_v = _mm_set1_ps(x[j]); __m128 xi_v = _mm_load_ps(&x[i]); __m128 rx_v = _mm_sub_ps(xj_v, xi_v); __m128 yj_v = _mm_set1_ps(y[j]); __m128 yi_v = _mm_load_ps(&y[i]); __m128 ry_v = _mm_sub_ps(yj_v, yi_v); __m128 zj_v = _mm_set1_ps(z[j]); __m128 zi_v = _mm_load_ps(&z[i]); __m128 rz_v = _mm_sub_ps(zj_v, zi_v); __m128 r2_v = _mm_mul_ps(rx_v, rx_v) + _mm_mul_ps(ry_v, ry_v) + _mm_mul_ps(rz_v, rz_v) + _mm_set1_ps(eps); __m128 r2inv_v = _mm_div_ps(_mm_set1_ps(1.0f),_mm_sqrt_ps(r2_v)); __m128 r6inv_1v = _mm_mul_ps(r2inv_v, r2inv_v); __m128 r6inv_v = _mm_mul_ps(r6inv_1v, r2inv_v); __m128 mj_v = _mm_set1_ps(m[j]); __m128 s_v = _mm_mul_ps(mj_v, r6inv_v); __m128 axi_v = _mm_load_ps(&ax[i]); __m128 ayi_v = _mm_load_ps(&ay[i]); __m128 azi_v = _mm_load_ps(&az[i]); __m128 srx_v = _mm_mul_ps(s_v, rx_v); __m128 sry_v = _mm_mul_ps(s_v, ry_v); __m128 srz_v = _mm_mul_ps(s_v, rz_v); axi_v = _mm_add_ps(axi_v, srx_v); ayi_v = _mm_add_ps(ayi_v, srx_v); azi_v = _mm_add_ps(azi_v, srx_v); _mm_store_ps(ax, axi_v); _mm_store_ps(ay, ayi_v); _mm_store_ps(az, azi_v); } } 

Je pense que l’idée principale est correcte mais il y a une / des erreur (s) quelque part car la réponse obtenue est incorrecte.

Merci d’avance.

Je pense que vos seuls bogues sont de simples fautes de frappe, pas des erreurs de logique, voir ci-dessous.

Ne pouvez-vous pas simplement utiliser la clang automatique de clang ? Ou avez-vous besoin d’utiliser gcc pour ce code? La vectorisation automatique vous permet de créer des versions SSE, AVX et (à l’avenir) AVX512 à partir de la même source sans modification. Les éléments insortingnsèques ne sont malheureusement pas adaptables à différentes tailles de vecteurs.

Sur la base de vos débuts en vectorisation, j’ai créé une version optimisée. Vous devriez l’essayer. Je suis curieux de savoir si c’est plus rapide que votre version comportant des corrections de bugs ou la version vectorisée automatique de Clang. 🙂


Cela semble faux:

 _mm_store_ps(ax, axi_v); _mm_store_ps(ay, ayi_v); _mm_store_ps(az, azi_v); 

Vous avez chargé depuis ax[i] , mais vous enregistrez maintenant dans ax[0] .


De plus, l’avertissement de variable inutilisée de clang a trouvé ce bogue:

 axi_v = _mm_add_ps(axi_v, srx_v); ayi_v = _mm_add_ps(ayi_v, srx_v); // should be sry_v azi_v = _mm_add_ps(azi_v, srx_v); // should be srz_v 

Comme je l’ai dit dans ma réponse à votre question précédente, vous devriez probablement échanger les boucles, donc le même axe [i + 0..3], ay [i + 0..3] et az [i + 0..3 ] sont utilisés en évitant ce chargement / stockage.

De plus, si vous n’utilisez pas rsqrtps + Newton-Raphson , vous devez utiliser la transformation que j’ai indiquée dans ma dernière réponse: divisez m[j] par sqrt(k2) ^3 . Il est inutile de diviser 1.0f par quelque chose utilisant un divps , puis de ne multiplier qu’une seule fois par la suite.

rsqrt pourrait ne pas être réellement gagnant, car le débit uop total pourrait être davantage un goulot d’étranglement que le débit ou la latence div / sqrt. trois multiplications + FMA + rsqrtps est significativement plus uops que sqrtps + divps. rsqrt est plus utile avec les vecteurs AVX 256b, car l’unité divide / sqrt n’est pas pleine largeur sur SnB vers Broadwell. Skylake a 12c de latence sqrtps ymm , comme pour xmm, mais le débit est toujours meilleur pour xmm (un par 3c au lieu de un par 6c).

clang et gcc utilisaient tous les deux rsqrtps / rsqrtss lors de la compilation de votre code avec -ffast-math . (seulement en utilisant la version packagée, bien sûr.)


Si vous n’échangez pas les boucles, vous devez lever manuellement tout ce qui ne dépend que de j hors de la boucle interne. Les compilateurs sont généralement bons à cela, mais il semble toujours judicieux de faire en sorte que la source reflète ce que vous attendez du compilateur. Cela aide à “voir” ce que la boucle est en train de faire.


Voici une version avec quelques optimisations par rapport à votre original:

  • Pour que gcc / clang fusionne le mul / ajoute dans FMA, j’ai utilisé -ffp-contract=fast . Cela obtient des instructions FMA pour un débit élevé sans utiliser -ffast-math . (Il y a beaucoup de parallélisme avec les trois accumulateurs séparés, donc la latence accrue de FMA par rapport aux addps ne devrait pas faire de mal. Je suppose que le débit en port0 / 1 est le facteur limitant ici.) Je pensais que gcc le faisait automatiquement , mais il semble que ce ne soit pas le cas ici sans -ffast-math .

  • Notez que v 3/2 = sqrt (v) 3 = sqrt (v) * v. Cela a une latence plus basse et moins d’instructions.

  • Interchangez les boucles et utilisez des charges de diffusion dans la boucle interne pour améliorer la localité (réduisez la bande passante requirejse de 4 ou 8 avec AVX). Chaque itération de la boucle interne ne lit que 4B de nouvelles données provenant de chacun des quatre tableaux source. (x, y, z et m). Donc, il fait beaucoup usage de chaque ligne de cache quand il fait chaud.

    L’utilisation de charges de diffusion dans la boucle interne signifie également que nous ax[i + 0..3] en parallèle, évitant ainsi la nécessité d’une sum horizontale, ce qui prend du code supplémentaire . (Voir une version précédente de cette réponse pour le code avec les boucles échangées, mais celui qui utilise des charges vectorielles dans la boucle interne, avec ssortingde = 16B .)

Il comstack bien pour Haswell avec gcc, en utilisant FMA . (Toujours en 128b, contrairement à la version 256b auto-vectorisée de Clang). La boucle interne ne contient que 20 instructions, dont 13 seulement sont des instructions FPU ALU nécessitant le port 0/1 de la famille Intel SnB. Il fait un code décent même avec SSE2 de base: pas de FMA, et a besoin de shufps pour les charges de diffusion, mais celles-ci ne se font pas concurrence pour les unités d’exécution avec add / mul.

 #include  void ffunc_sse128(float *ressortingct ax, float *ressortingct ay, float *ressortingct az, const float *x, const float *y, const float *z, int N, float eps, const float *ressortingct m) { for (int i = 0; i < N; i+=4) { __m128 xi_v = _mm_load_ps(&x[i]); __m128 yi_v = _mm_load_ps(&y[i]); __m128 zi_v = _mm_load_ps(&z[i]); // vector accumulators for ax[i + 0..3] etc. __m128 axi_v = _mm_setzero_ps(); __m128 ayi_v = _mm_setzero_ps(); __m128 azi_v = _mm_setzero_ps(); // AVX broadcast-loads are as cheap as normal loads, // and data-reuse meant that stand-alone load instructions were used anyway. // so we're not even losing out on folding loads into other insns // An inner-loop stride of only 4B is a huge win if memory / cache bandwidth is a bottleneck // even without AVX, the shufps instructions are cheap, // and don't compete with add/mul for execution units on Intel for (int j = 0; j < N; j++) { __m128 xj_v = _mm_set1_ps(x[j]); __m128 rx_v = _mm_sub_ps(xj_v, xi_v); __m128 yj_v = _mm_set1_ps(y[j]); __m128 ry_v = _mm_sub_ps(yj_v, yi_v); __m128 zj_v = _mm_set1_ps(z[j]); __m128 rz_v = _mm_sub_ps(zj_v, zi_v); __m128 mj_v = _mm_set1_ps(m[j]); // do the load early // sum of squared differences __m128 r2_v = _mm_set1_ps(eps) + rx_v*rx_v + ry_v*ry_v + rz_v*rz_v; // GNU extension /* __m128 r2_v = _mm_add_ps(_mm_set1_ps(eps), _mm_mul_ps(rx_v, rx_v)); r2_v = _mm_add_ps(r2_v, _mm_mul_ps(ry_v, ry_v)); r2_v = _mm_add_ps(r2_v, _mm_mul_ps(rz_v, rz_v)); */ // rsqrt and a Newton-Raphson iteration might have lower latency // but there's enough surrounding instructions and cross-iteration parallelism // that the single-uop sqrtps and divps instructions prob. aren't be a bottleneck __m128 r2sqrt = _mm_sqrt_ps(r2_v); __m128 r6sqrt = _mm_mul_ps(r2_v, r2sqrt); // v^(3/2) = sqrt(v)^3 = sqrt(v)*v __m128 s_v = _mm_div_ps(mj_v, r6sqrt); __m128 srx_v = _mm_mul_ps(s_v, rx_v); __m128 sry_v = _mm_mul_ps(s_v, ry_v); __m128 srz_v = _mm_mul_ps(s_v, rz_v); axi_v = _mm_add_ps(axi_v, srx_v); ayi_v = _mm_add_ps(ayi_v, sry_v); azi_v = _mm_add_ps(azi_v, srz_v); } _mm_store_ps(&ax[i], axi_v); _mm_store_ps(&ay[i], ayi_v); _mm_store_ps(&az[i], azi_v); } } 

J'ai aussi essayé une version avec rcpps , mais IDK si ce sera plus rapide. Notez qu'avec -ffast-math , gcc et clang convertiront la division en rcpps + une itération de Newton. (Pour une raison quelconque, ils ne convertissent pas 1.0f/sqrtf(x) en rsqrt + Newton, même dans une fonction autonome). Clang fait un meilleur travail en utilisant FMA pour l'étape d'itération.

 #define USE_RSQRT #ifndef USE_RSQRT // even with -mrecip=vec-sqrt after -ffast-math, this still does sqrt(v)*v, then rcpps __m128 r2sqrt = _mm_sqrt_ps(r2_v); __m128 r6sqrt = _mm_mul_ps(r2_v, r2sqrt); // v^(3/2) = sqrt(v)^3 = sqrt(v)*v __m128 s_v = _mm_div_ps(mj_v, r6sqrt); #else __m128 r2isqrt = rsqrt_float4_single(r2_v); // can't use the sqrt(v)*v sortingck, unless we either do normal sqrt first then rcpps // or rsqrtps and rcpps. Maybe it's possible to do a Netwon Raphson iteration on that product // instead of refining them both separately? __m128 r6isqrt = r2isqrt * r2isqrt * r2isqrt; __m128 s_v = _mm_mul_ps(mj_v, r6isqrt); #endif