Le moyen le plus rapide de multiplier un tableau de int64_t?

Je souhaite vectoriser la multiplication de deux tableaux alignés sur la mémoire. Je n’ai trouvé aucun moyen de multiplier 64 * 64 bits dans AVX / AVX2; Y a-t-il un moyen plus rapide de faire cela?

Remarque: je ne veux pas enregistrer le résultat de moitié élevé de chaque multiplication.

void multiply_vex(long *Gi_vec, long q, long *Gj_vec){ int i; __m256i data_j, data_i; __uint64_t *ptr_J = (__uint64_t*)&data_j; __uint64_t *ptr_I = (__uint64_t*)&data_i; for (i=0; i<BASE_VEX_STOP; i+=4) { data_i = _mm256_load_si256((__m256i*)&Gi_vec[i]); data_j = _mm256_load_si256((__m256i*)&Gj_vec[i]); ptr_I[0] -= ptr_J[0] * q; ptr_I[1] -= ptr_J[1] * q; ptr_I[2] -= ptr_J[2] * q; ptr_I[3] -= ptr_J[3] * q; _mm256_store_si256((__m256i*)&Gi_vec[i], data_i); } for (; i<BASE_DIMENSION; i++) Gi_vec[i] -= Gj_vec[i] * q; } 

UPDATE: J’utilise la microarchitecture Haswell avec les deux compilateurs ICC / GCC. Donc, les deux AVX et AVX2 va bien. Je remplace le -= par le C insortingsique _mm256_sub_epi64 après la multiplication en boucle-dérouler, où il obtient une certaine accélération. Actuellement, c’est ptr_J[0] *= q; ... ptr_J[0] *= q; ...

J’utilise __uint64_t mais c’est une erreur . Le bon type de données est __int64_t .

    Vous semblez supposer que long est 64 bits dans votre code, mais utilisez ensuite __uint64_t également. En 32 bits, l’ ABI x32 et, sous Windows, long est un type 32 bits. Votre titre mentionne long long , mais ensuite votre code l’ignore. Je me demandais depuis un moment si votre code supposait que 32 bits était long .

    Vous vous tirez complètement dans le pied en utilisant des charges AVX256, puis en aliasant un pointeur sur le __m256i pour effectuer des opérations scalaires. gcc abandonne juste et vous donne le code terrible que vous avez demandé: charge de vecteur puis un tas d’ extract et d’ insert instructions. Votre façon d’écrire cela signifie que les deux vecteurs doivent également être décompressés pour vpsubq sous- vpsubq scalaire, au lieu d’utiliser vpsubq .

    Les processeurs x86 modernes disposent d’un cache L1 très rapide pouvant gérer deux opérations par horloge . (Haswell et plus tard: deux charges et un magasin par horloge). Effectuer plusieurs charges scalaires à partir de la même ligne de cache est préférable à un chargement et un décompactage vectoriels. (Cependant, la planification imparfaite du fonctionnement optimisé réduit le débit à environ 84%: voir ci-dessous)


    gcc 5.3 -O3 -march = haswell (explorateur du compilateur Godbolt) auto-vectorise assez bien une implémentation scalaire simple. Lorsque AVX2 n’est pas disponible, gcc continue de s’auto-vectoriser avec des vecteurs de 128 bits: sous Haswell, la vitesse du code 64 bits scalaire idéal sera en fait environ la moitié. (Voir l’parsing perf ci-dessous, mais remplacez 2 éléments par vecteur au lieu de 4).

     #include  // why not use this like a normal person? #define BASE_VEX_STOP 1024 #define BASE_DIMENSION 1028 // ressortingct lets the comstackr know the arrays don't overlap, // so it doesn't have to generate a scalar fallback case void multiply_simple(uint64_t *ressortingct Gi_vec, uint64_t q, const uint64_t *ressortingct Gj_vec){ for (intptr_t i=0; i 

    Boucle intérieure:

     .L4: vmovdqu ymm1, YMMWORD PTR [r9+rax] # MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B] add rcx, 1 # ivtmp.30, vpsrlq ymm0, ymm1, 32 # tmp174, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], vpmuludq ymm2, ymm1, ymm3 # tmp173, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], vect_cst_.25 vpmuludq ymm0, ymm0, ymm3 # tmp176, tmp174, vect_cst_.25 vpmuludq ymm1, ymm4, ymm1 # tmp177, tmp185, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B] vpaddq ymm0, ymm0, ymm1 # tmp176, tmp176, tmp177 vmovdqa ymm1, YMMWORD PTR [r8+rax] # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B] vpsllq ymm0, ymm0, 32 # tmp176, tmp176, vpaddq ymm0, ymm2, ymm0 # vect__13.24, tmp173, tmp176 vpsubq ymm0, ymm1, ymm0 # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24 vmovdqa YMMWORD PTR [r8+rax], ymm0 # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26 add rax, 32 # ivtmp.32, cmp rcx, r10 # ivtmp.30, bnd.14 jb .L4 #, 

    Si vous le souhaitez, remettez-le à la source, mais il sera beaucoup plus facile de laisser le compilateur s'auto-vectoriser. Je n'ai pas essayé de l'parsingr pour voir si c'était optimal.

    Si vous ne comstackz généralement pas avec -O3 , vous pouvez utiliser #pragma omp simd avant la boucle (et -fopenmp ).

    Bien sûr, au lieu d’un épilogue scalaire, c’est prob. accélérez le chargement des 32 derniers bits de Gj_vec et stockez-les dans les derniers 32B de Gi_vec, en faisant éventuellement chevauchement avec le dernier magasin de la boucle. (Un repli scalaire est toujours nécessaire si les tableaux sont plus petits que 32B.)


    Version insortingnsèque de vecteur améliorée pour Haswell

    De mes commentaires sur la réponse de Z Boson. Basé sur le code de la bibliothèque de classes vectorielles d'Agner Fog .

    La version d'Agner Fog enregistre une instruction mais des goulots d'étranglement sur le port de lecture aléatoire en utilisant phadd + pshufd, où j'utilise psrlq / paddq / pand.

    Etant donné que l’un de vos opérandes est constant, assurez-vous de passer set1(q) sous la forme b , et non a , afin que le shuffle "bswap" puisse être levé.

     // replace hadd -> shuffle (4 uops) with shift/add/and (3 uops) // The constant takes 2 insns to generate outside a loop. __m256i mul64_haswell (__m256i a, __m256i b) { // instruction does not exist. Split into 32-bit multiplies __m256i bswap = _mm256_shuffle_epi32(b,0xB1); // swap H<->L __m256i prodlh = _mm256_mullo_epi32(a,bswap); // 32 bit L*H products // or use pshufb instead of psrlq to reduce port0 pressure on Haswell __m256i prodlh2 = _mm256_srli_epi64(prodlh, 32); // 0 , a0Hb0L, 0, a1Hb1L __m256i prodlh3 = _mm256_add_epi32(prodlh2, prodlh); // xxx, a0Lb0H+a0Hb0L, xxx, a1Lb1H+a1Hb1L __m256i prodlh4 = _mm256_and_si256(prodlh3, _mm256_set1_epi64x(0x00000000FFFFFFFF)); // zero high halves __m256i prodll = _mm256_mul_epu32(a,b); // a0Lb0L,a1Lb1L, 64 bit unsigned products __m256i prod = _mm256_add_epi64(prodll,prodlh4); // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32 return prod; } 

    Voir sur Godbolt .

    Notez que cela n'inclut pas la soustraction finale, mais seulement la multiplication.

    Cette version devrait fonctionner un peu mieux sur Haswell que la version autovectorisée de gcc. (comme peut-être un vecteur par 4 cycles au lieu d'un vecteur par 5 cycles, goulot d'étranglement sur le débit du port0. Je n'ai pas pris en compte d'autres goulots d'étranglement pour le problème complet, car c'était un ajout tardif à la réponse.)

    Une version AVX1 (deux éléments par vecteur) serait nulle, et serait probablement encore pire que le scalaire 64 bits. Ne le faites pas sauf si vous avez déjà vos données dans des vecteurs et que vous voulez obtenir le résultat dans un vecteur (extraire en scalaire et vice versa pourrait ne pas en valoir la peine).


    Analyse approfondie du code autovectorisé de GCC (pas la version insortingnsèque)

    Contexte: voir les tableaux insn et le guide microarch d'Agner Fog , ainsi que d'autres liens dans le wiki des balises x86

    Jusqu'en AVX512 (voir ci-dessous), ce n'est probablement qu'un peu plus rapide que le code scalaire 64 bits: imul r64, m64 a un débit d'un par horloge sur les processeurs Intel (mais un sur quatre par la famille AMD Bulldozer). load / imul / sub-with-memory-dest correspond à 4 uops de domaine fondu sur des processeurs Intel (avec un mode d'adressage pouvant effectuer une micro-fusion, gcc ne pouvant pas être utilisé). La largeur du pipeline est de 4 UOP de domaine fondu par horloge. Par conséquent, même un grand déroulé ne peut pas le faire émettre à un par horloge. Avec suffisamment de déroulage, nous allons créer un goulot d'étranglement sur le débit de chargement / stockage. Haswell permet deux chargements et un magasin par horloge, mais , selon le manuel d’Intel, le débit des ports de chargement volés d’adresses par adresses abaissera le débit à environ 81/96 = 84% .

    Alors peut-être que le meilleur moyen pour Haswell serait de charger et de multiplier avec scalaire (2 uops), puis vmovq / pinsrq / vinserti128 afin que vous puissiez faire la soustraction avec un vpsubq . C'est 8 uops pour charger et multiplier les 4 scalaires, 7 shops pour obtenir les données dans un __m256i (2 (movq) + 4 (pinsrq est 2 uops) + 1 vinserti128), et 3 autres uops pour faire un vecteur load / vpsubq / vector le magasin. Il s’agit donc de 18 Uops à domaine fondu pour 4 multiplications (4,5 cycles à émettre), mais 7 Uops shuffle (7 cycles à exécuter). Donc nvm, ce n'est pas bon comparé au scalaire pur.


    Le code autovectorisé utilise 8 instructions ALU de vecteur pour chaque vecteur de quatre valeurs. Sur Haswell, 5 de ces uops (multiplications et décalages) ne peuvent s'exécuter que sur le port 0. Ainsi, quel que soit le déroulement de cet algorithme, il obtiendra au mieux un vecteur pour 5 cycles (c'est-à-dire une multiplication par 5/4 cycles).

    Les décalages pourraient être remplacés par pshufb (port 5) pour déplacer les données et les décalés par zéros. (Les autres shuffles ne prennent pas en charge la réduction à zéro au lieu de copier un octet depuis l'entrée et il n'y a pas de zéros connus dans l'entrée que nous pourrions copier.)

    paddq / psubq peut fonctionner sur les ports 1/5 sur Haswell ou p015 sur Skylake.

    Skylake exécute pmuludq et les décalages vectoriels de comptage immédiat sur p01, de sorte qu’il pourrait en théorie gérer un débit d’un vecteur par max (5/2, 8/3, 11/4) = 11/4 = 2,75 cycles. Cela engendre donc des goulots d'étranglement sur le débit total d'uop du domaine fondu (y compris les 2 charges vectorielles et 1 magasin de vecteurs). Donc, un peu de déroulement de boucle aidera. Il est probable que des conflits de ressources liés à une planification imparfaite entraîneront un goulot d'étranglement à un peu moins de 4 Uops de domaine fondu par horloge. La surcharge de la boucle peut éventuellement être exécutée sur le port 6, qui ne peut gérer que certaines opérations scalaires, y compris les opérations d’ add et de comparaison, laissant ainsi les ports 0/1/5 pour les opérations ALU vectorielles, car elles sont proches de la saturation (8/3 = 2,666 horloges). Les ports de chargement / stockage sont loin d’être saturés, cependant.

    Ainsi, Skylake peut théoriquement gérer un vecteur par 2,75 cycles (surcharge de la boucle), ou un multiplié par environ 0,7 cycles , par rapport à la meilleure option de Haswell (un par 1,2 cycles en théorie avec scalaire, ou un par 1,25 cycles en théorie avec des vecteurs ). Le cycle scalaire un pour environ 1,2 cycle nécessiterait probablement une boucle asm réglée manuellement, car les compilateurs ne savent pas comment utiliser un mode d'adressage à un registre pour les magasins et un mode d'adressage à deux registres pour les charges ( dst + (src-dst) et incrémentation dst ).

    De plus, si vos données ne sont pas à jour dans le cache L1, le travail effectué avec moins d'instructions permet à l'interface client de devancer les unités d'exécution et de commencer les chargements avant que les données ne soient nécessaires. La pré-extraction matérielle ne traverse pas les lignes de page, aussi une boucle vectorielle battra-t-elle probablement scalar en pratique pour les tableaux de grande taille, et peut-être même pour les tableaux plus petits .


    AVX-512DQ introduit un vecteur multiplier 64bx64b-> 64b

    Si vous ajoutez -mavx512dq gcc peut se vectoriser -mavx512dq .

     .L4: vmovdqu64 zmm0, ZMMWORD PTR [r8+rax] # vect__11.23, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B] add rcx, 1 # ivtmp.30, vpmullq zmm1, zmm0, zmm2 # vect__13.24, vect__11.23, vect_cst_.25 vmovdqa64 zmm0, ZMMWORD PTR [r9+rax] # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B] vpsubq zmm0, zmm0, zmm1 # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24 vmovdqa64 ZMMWORD PTR [r9+rax], zmm0 # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26 add rax, 64 # ivtmp.32, cmp rcx, r10 # ivtmp.30, bnd.14 jb .L4 #, 

    Ainsi, AVX512DQ (qui devrait faire partie de Skylake multi-socket Xeon (Purley) vers 2017 ) donnera une accélération beaucoup plus grande que 2x (à partir de vecteurs plus larges) si ces instructions sont traitées en pipeline une par horloge.

    Mise à jour: Skylake-AVX512 (ou SKL-X ou SKL-SP) exécute VPMULLQ à un cycle et demi pour les vecteurs xmm, ymm ou zmm. C'est 3 uops avec 15c de latence. (Avec peut-être une latence supplémentaire de 1c pour la version zmm, si ce n'est pas un problème de mesure dans les résultats AIDA .)

    vpmullq est beaucoup plus rapide que tout ce que vous pouvez créer à partir de vpmullq 32 bits, il est donc très utile d'avoir une instruction à ce sujet, même si les processeurs actuels ne disposent pas de matériel à multiplication de vecteurs de 64 bits. (Ils utilisent vraisemblablement les multiplicateurs de mantisse dans les unités FMA.)

    Si les opérations SIMD 64bx64b à 64b (moins) vous intéressent, voici les solutions AVX et AVX2 de la bibliothèque de classes vectorielles d’Agner Fog. Je voudrais les tester avec des tableaux et voir comment il se compare à ce que fait GCC avec une boucle générique telle que celle de la réponse de Peter Cordes.

    AVX (utilisez SSE – vous pouvez toujours comstackr avec -mavx pour obtenir le codage -mavx ).

     // vector operator * : multiply element by element static inline Vec2q operator * (Vec2q const & a, Vec2q const & b) { #if INSTRSET >= 5 // SSE4.1 supported // instruction does not exist. Split into 32-bit multiplies __m128i bswap = _mm_shuffle_epi32(b,0xB1); // b0H,b0L,b1H,b1L (swap H<->L) __m128i prodlh = _mm_mullo_epi32(a,bswap); // a0Lb0H,a0Hb0L,a1Lb1H,a1Hb1L, 32 bit L*H products __m128i zero = _mm_setzero_si128(); // 0 __m128i prodlh2 = _mm_hadd_epi32(prodlh,zero); // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0 __m128i prodlh3 = _mm_shuffle_epi32(prodlh2,0x73); // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L __m128i prodll = _mm_mul_epu32(a,b); // a0Lb0L,a1Lb1L, 64 bit unsigned products __m128i prod = _mm_add_epi64(prodll,prodlh3); // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32 return prod; #else // SSE2 int64_t aa[2], bb[2]; a.store(aa); // split into elements b.store(bb); return Vec2q(aa[0]*bb[0], aa[1]*bb[1]); // multiply elements separetely #endif } 

    AVX2

     // vector operator * : multiply element by element static inline Vec4q operator * (Vec4q const & a, Vec4q const & b) { // instruction does not exist. Split into 32-bit multiplies __m256i bswap = _mm256_shuffle_epi32(b,0xB1); // swap H<->L __m256i prodlh = _mm256_mullo_epi32(a,bswap); // 32 bit L*H products __m256i zero = _mm256_setzero_si256(); // 0 __m256i prodlh2 = _mm256_hadd_epi32(prodlh,zero); // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0 __m256i prodlh3 = _mm256_shuffle_epi32(prodlh2,0x73); // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L __m256i prodll = _mm256_mul_epu32(a,b); // a0Lb0L,a1Lb1L, 64 bit unsigned products __m256i prod = _mm256_add_epi64(prodll,prodlh3); // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32 return prod; } 

    Ces fonctions fonctionnent pour les entiers 64 bits signés et non signés. Dans votre cas, puisque q est constant dans la boucle, vous n'avez pas besoin de recalculer certaines choses à chaque itération, mais votre compilateur le comprendra probablement de toute façon.