Comment initialiser un vecteur SIMD avec une plage de 0 à N?

J’ai la fonction suivante pour laquelle j’essaye d’écrire une version AXV:

void hashids_shuffle(char *str, size_t str_length, char *salt, size_t salt_length) { size_t i, j, v, p; char temp; if (!salt_length) { return; } for (i = str_length - 1, v = 0, p = 0; i > 0; --i, ++v) { v %= salt_length; p += salt[v]; j = (salt[v] + v + p) % i; temp = str[i]; str[i] = str[j]; str[j] = temp; } } 

J’essaie de vectoriser v %= salt_length; .
Je souhaite initialiser un vecteur contenant des nombres compris entre 0 et str_length afin d’utiliser _mm_rem_epu64 de SVML afin de calculer v pour chaque itération de boucle.
Comment initialiser le vecteur correctement?

Demander comment initialiser un vecteur revient en gros à demander un tutoriel. Google place certains des guides d’Intel sur l’utilisation d’insortingnsèques. Je vais faire comme si la question n’était pas anodine, et répondre à propos d’essayer de mettre en œuvre cette fonction efficacement . Ce n’est certainement pas le genre de fonction que vous essayez de vectoriser en tant que débutant total.

Voir le wiki des balises x86 pour des liens vers des documents, en particulier. Guide insortingnsèque d’Intel. Voir le wiki de la balise sse pour un lien vers une très bonne introduction à la programmation SIMD avec les éléments insortingnsèques SSE, et comment utiliser efficacement SIMD , entre autres liens.

Sommaire sommaire :

  • Unrolling / vectorizing prend en charge gratuitement v % salt_length .
  • Comment vous pourriez vectoriser v++; v %= loop_invariant; v++; v %= loop_invariant; si ce n’était pas une puissance de 2 ou une constante de compilation. Inclut une réponse à la question sur l’utilisation de _mm_set_epi8 ou d’autres moyens d’initialiser un vecteur à cette fin.
  • Comment commencer à vectoriser une boucle complexe comme ceci: commencez par dérouler pour trouver des dépendances en série.
  • une version non testée de la fonction complète qui vectorise tout sauf le % i et le swap. (c’est-à-dire que toutes les opérations qui étaient bon marché, comme vous l’avez demandé).

    • (v + salt[v] + p) (et tout ce qui le vpaddw ) vectorise deux instructions vpaddw . La configuration préfixe-sum en dehors de la boucle pour la vectorisation de p était délicate, mais je l’ai également vectorisée.

    • La grande majorité du temps d’exécution de la fonction sera dans la boucle interne scalaire recouvrant un vecteur d’éléments j , goulot d’étranglement sur div (ou quoi que SVML puisse faire) et / ou des erreurs de cache avec de très grandes chaînes.


La totalité de la boucle ne peut pas être vectorisée facilement, car les échanges avec des indices pseudo-aléatoires créent une dépendance en série imprévisible. Utiliser AVX512 pour rassembler + mélanger + disperser, avec AVX512CD pour rechercher des masques de bits en conflit, est peut-être possible, mais il faudrait que ce soit une question distincte. Je ne sais pas à quel point il serait difficile de faire cela efficacement, ou si vous finissiez souvent par répéter un mélange aléatoire de vecteur, ne réalisant des progrès que dans un élément non conflictuel.


Puisque salt_length = sizeof(size_t) est une constante de compilation et une puissance inférieure de 2 à la longueur de votre vecteur, v++ et v%=salt_length ne nécessitent aucun code dans la boucle , et se produisent gratuitement en tant que effet de dérouler efficacement la boucle pour faire plusieurs valeurs v en parallèle.

(L’utilisation d’une taille de sel dépendant de la plate-forme signifie qu’une version 32 bits ne pourra pas traiter les données créées avec du sel 64 bits. Même l’ABI x32 a une size_t 32 bits, donc passer à uint64_t semblerait logique, sauf si vous n’avez jamais besoin de partager des hachures salées entre des machines.

Dans la boucle scalaire, v suit un motif répétitif de 0 1 2 3 0 1 2 3 … (ou 0..7 pour 64 bits). En code vectoriel, nous faisons peut-être 8 valeurs v à la fois avec 4 éléments dans un vecteur de 32 octets, ou 16 itérations à la fois avec 2 éléments.

Donc, v devient un vecteur constant invariant des boucles. Il est intéressant de noter que salt[v] , nous n’avons donc jamais besoin de faire des recherches dans la table des sels dans la boucle. En fait, v+salt[v] peut être pré-calculé pour scalaire et vecteur.

La version scalaire doit pré-calculer v+salt[v] et se dérouler par 4 ou 8 également, en supprimant la recherche LUT afin que tout le débit de mémoire / cache soit disponible pour les échanges réels. Le compilateur ne le fera probablement pas pour vous. Vous devrez donc probablement vous dérouler manuellement et écrire du code supplémentaire pour gérer le dernier nombre impair d’octets de chaîne. (Sans dérouler, vous pouvez toujours pré-calculer une table de recherche de v+salt[v] , avec un type suffisamment large pour ne pas envelopper).

Même s’assurer que salt_length est connu au moment de la compilation permettrait également un code bien meilleur. v %= comstack_time_constant est moins cher qu’un div insn, et extrêmement bon marché lorsqu’il est une puissance de 2. (Il se transforme simplement en v &= 7 ). Le compilateur peut éventuellement le faire pour vous si la version scalaire peut être intégrée, ou si vous avez utilisé salt_length = sizeof(size_t) au lieu de le transmettre sous la forme d’une argument de fonction.


Si vous ne connaissiez pas déjà salt_length: c’est -à- dire ce que @harold suggérait avant de révéler les informations critiques sur salt_length:

Puisque nous connaissons v < salt_length pour commencer, nous n’avons besoin que d’une seule v -= salt_length pour le replacer dans la bonne plage et maintenir cet invariant. C'est ce qu'on appelle une opération de "réduction de force", car la soustraction est une opération plus faible (et moins chère) que la division.

 // The scalar loop would benefit from this transformation, too. // or better, unroll the scalar loop by 8 so everything becomes constant v++; if( v >= salt_length) v-= salt_length; 

Pour vectoriser simplement ceci: supposons que tout ce que nous soaps est salt_length <= 16 , nous pouvons donc utiliser un vecteur de 32 valeurs uint8_t. (Et nous pouvons utiliser pshufb pour vectoriser la recherche salt[v] LUT).

 // untested // Vectorizing v++; v %= unknown_loop_invariant_value; if (!salt_length) return; assert(salt_length <= 16); // so we can use pshufb for the salt[v] step __m256i vvec = _mm256_setr_epi8( // setr: lowest element first, unlike set 0%salt_length, 1%salt_length, 2%salt_length, 3%salt_length, 4%salt_length, 5%salt_length, 6%salt_length, 7%salt_length, 8%salt_length, 9%salt_length, 10%salt_length, 11%salt_length, 12%salt_length, 13%salt_length, 14%salt_length, 15%salt_length, 16%salt_length, 17%salt_length, 18%salt_length, 19%salt_length, 20%salt_length, 21%salt_length, 22%salt_length, 23%salt_length, 24%salt_length, 25%salt_length, 26%salt_length, 27%salt_length, 28%salt_length, 29%salt_length, 30%salt_length, 31%salt_length); __m256i v_increment = _mm256_set1_epi8(32 % salt_length); __m256i vmodulus = _mm256_set1_epi8(salt_length); // salt_lut = _mm256_set1_epi64x(salt_byval); // for known salt length. (pass it by val in a size_t arg, instead of by char*). // duplicate the salt into both lanes of a vector. Garbage beyond salt_length isn't looked at. __m256i salt_lut = _mm256_broadcastsi128_si256(_mm_loadu_si128(salt)); // nevermind that this could segfault if salt is short and at the end of a page. //__m256i v_plus_salt_lut = _mm256_add_epi8(vvec, salt_lut); // not safe with 8-bit elements: could wrap // We could use 16-bit elements and AVX512 vpermw (or vpermi2w to support longer salts) for (...) { vvec = _mm256_add_epi8(vvec, v_increment); // ++v; // if(!(salt_length > v)) { v-= salt_length; } __m256i notlessequal = _mm256_cmpgt_epi8(vmodulus, vvec); // all-ones where salt_length > v. // all-zero where salt_length <= v, where we need to subtract salt_length __m256i conditional_sub = _mm256_and_si256(notlessequal, vmodulus) vvec = _mm256_sub_epi8(vvec, conditional_sub); // subtract 0 or salt_length // salt[v] lookup: __m256i saltv = _mm256_shuffle_epi8(salt_lut, vvec); // salt[v] // then maybe pmovzx and vextracti128+pmovzx to zero-extend to 16-bit elements? Maybe vvec should only be a 16-bit vector? // or unpack lo/hi with zeros (but that behaves differently from pmovzx at the lane boundary) // or have vvec already holding 16-bit elements with the upper half of each one always zero. mask after the pshufb to re-zero, // or do something clever with `vvec`, `v_increment` and `vmodulus` so `vvec` can have `0xff` in the odd bytes, so pshufb zeros those elements. } 

Bien sûr, si nous savions que salt_length était une puissance de 2, nous aurions simplement masqué tous les bits sauf les bits bas pertinents dans chaque élément:

 vvec = _mm256_add_epi8(vvec, _mm256_set1_epi8(salt_length)); // ++v; vvec = _mm256_and_si256(vvec, _mm256_set1_epi8(salt_length - 1)); // v &= salt_length - 1; // aka v%=salt_length; 

Lorsque nous avons commencé avec la mauvaise taille d’élément, nous nous sums rendu compte que ne vectoriser qu’une ligne à la fois était une mauvaise idée, car nous devons maintenant modifier tout le code que nous avons déjà écrit pour utiliser des éléments plus larges, nécessitant parfois une stratégie différente. faire la même chose.

Bien sûr, vous devez commencer par un aperçu, mental ou écrit, de la manière dont vous pourriez effectuer chaque étape séparément. C'est en train de réfléchir à cela que vous voyez comment les différentes parties pourraient s'emboîter.

Pour les boucles complexes, une première étape utile consiste peut-être à dérouler manuellement la boucle scalaire. Cela aidera à trouver des dépendances de série et des choses qui simplifient avec le déroulement.


(stuff) % i : La partie difficile

Nous avons besoin d'éléments suffisamment larges pour contenir la valeur maximale de i , car i n'est pas une puissance de 2, ni constante, de sorte que l'opération modulo demande du travail. Toute largeur est un gaspillage et réduit notre débit. Si nous pouvions vectoriser tout le rest de la boucle, il pourrait même être intéressant de spécialiser la fonction avec différentes versions pour différentes plages de str_length . (Ou peut-être en boucle avec des éléments 64b jusqu'à i <= UINT32_MAX, puis en boucle jusqu'à i <= UINT16_MAX, etc.). Si vous savez que vous n'avez pas besoin de gérer les chaînes> 4GiB, vous pouvez accélérer le cas courant en utilisant uniquement des mathématiques 32 bits. (La division 64 bits est plus lente que la division 32 bits, même lorsque les bits supérieurs sont tous nuls).

En fait, je pense que nous avons besoin d’éléments aussi larges que le maximum p , puisqu’il continue à s’accumuler pour toujours (jusqu’à ce que le code scalaire initial se termine à 2 ^ 64). Contrairement à un module constant, nous ne pouvons pas simplement utiliser p%=i pour le contrôler, même si modulo est dissortingbutif. (123 % 33) % (33-16) != 123 % (33-16) . Même s’aligner sur 16 n’aide pas: 12345% 32! = 12345% 48% 32

Cela rendra rapidement p trop grand pour une soustraction conditionnelle répétée de i (jusqu'à ce que le masque de condition soit entièrement faux), même pour des valeurs relativement grandes de i .

Il existe des astuces pour modulo avec des constantes entières connues (voir http://libdivide.com/ ), mais AFAIK qui travaille sur l’inverse modulaire multiplicatif pour un diviseur proche (même avec une foulée de puissance double comme 16) n’est plus facile que pour un numéro totalement séparé. Nous ne pouvions donc pas ajuster à bas prix les constantes pour le prochain vecteur de valeurs i .

La loi des petits nombres rend peut-être utile de décoller les dernières itérations de vecteurs, avec des vecteurs pré-calculés d’inverses modulaires multiplicatifs afin que le % i puisse être effectué avec des vecteurs. Une fois que nous approchons de la fin de la chaîne, il fait probablement chaud dans le cache L1, nous sums donc totalement goulot d’étranglement sur div , pas sur les chargements / magasins d’échange. Pour cela, nous utiliserions peut-être un prolog pour atteindre une valeur i multiple de 16, de sorte que les derniers vecteurs à l'approche de i = 0 ont toujours le même alignement des valeurs i . Sinon, nous aurions une LUT de constantes pour une plage de valeurs i et ferions simplement des charges non alignées à partir de celle-ci. Cela signifie que nous n'avons pas besoin de faire pivoter salt_v et p .

Une conversion en FP serait peut-être utile, car les processeurs Intel récents (en particulier Skylake) possèdent un matériel de division FP très puissant avec un pipeline important (débit: rapport de latence). Si nous pouvons obtenir des résultats exacts avec le bon choix d'arrondi, ce serait formidable. ( float et double peuvent représenter exactement n'importe quel entier jusqu'à environ la taille de leur mantisse.)

Je suppose que cela vaut la peine d'essayer _mm_rem_epu16 d'Intel (avec un vecteur de i que vous décrémentez avec un vecteur de set1(16) ). S'ils utilisent la PF pour obtenir des résultats précis, c'est génial. Si elle décompresse simplement scalar et effectue une division entière, cela perdrait du temps à récupérer des valeurs dans un vecteur.

Quoi qu'il en soit, la solution la plus simple consiste certainement à itérer des éléments vectoriels avec une boucle scalaire. Tant que vous n'utiliserez pas AVX512CD pour les échanges, cela semble raisonnable, mais il s'agit probablement d'un ordre de grandeur plus lent que les échanges, s'ils sont tous actifs dans le cache N1.


(Untest) version partiellement vectorisée de la fonction:

Voici le code de l'explorateur du compilateur Godbolt , avec ses commentaires de conception complets, y compris les schémas que j'ai réalisés lors de la définition de l'algo préfixe-sum SIMD. Je me suis finalement souvenu que j'avais vu une version plus étroite de ceci comme bloc de construction dans la réponse à la sum des préfixes SSE en virgule flottante de @ ZBoson , mais ce n'est qu'après l'avoir essentiellement réinventée moi-même.

 // See the godbolt link for full design-notes comments // comments about what makes nice ASM or not. #include  #include  #include  #include  static inline __m256i init_p_and_increment(size_t salt_length, __m256i *p_increment, __m256i saltv_u16, __m128i saltv_u8) { // return initial p vector (for first 16 i values). // return increment vector by reference. if (salt_length == 4) { assert(0); // unimplemented // should be about the same as length == 8. Can maybe factor out some common parts, like up to psum2 } else { assert(salt_length == 8); // SIMD prefix sum for n elements in a vector in O(log2(n)) steps. __m128i sv = _mm256_castsi256_si128(saltv_u16); __m128i pshift1 = _mm_bslli_si128(sv, 2); // 1 elem (uint16_t) __m128i psum1 = _mm_add_epi16(pshift1, sv); __m128i pshift2 = _mm_bslli_si128(psum1, 4); // 2 elem __m128i psum2 = _mm_add_epi16(pshift2, psum1); __m128i pshift3 = _mm_bslli_si128(psum2, 8); // 4 elem __m128i psum3 = _mm_add_epi16(pshift3, psum2); // p_initial low 128. 2^3 = 8 elements = salt_length // psum3 = the repeating pattern of p values. Later values just add sum(salt[0..7]) to every element __m128i p_init_low = psum3; __m128i sum8_low = _mm_sad_epu8(saltv_u8, _mm_setzero_si128()); // sum(s0..s7) in each 64-bit half // alternative: // sum8_low = _mm_bsrli_si128(p_init_low, 14); // has to wait for psum3 to be ready: lower ILP than doing psadbw separately __m256i sum8 = _mm256_broadcastw_epi16(sum8_low); *p_increment = _mm256_slli_epi16(sum8, 1); // set1_epi16(2*sum(salt[0..7])) __m128i p_init_high = _mm_add_epi16(p_init_low, _mm256_castsi256_si128(sum8)); __m256i p_init = _mm256_castsi128_si256(p_init_low); p_init = _mm256_inserti128_si256(p_init, p_init_high, 1); // not supported by gcc _mm256_set_m128i(p_init_high, psum3); return p_init; } } void hashids_shuffle_simd(char *ressortingct str, size_t str_length, size_t salt_byval) { //assert(salt_length <= 16); // so we can use pshufb for the salt[v] step for non-constant salt length. // platform-dependent salt size seems weird. Why not uint64_t? size_t salt_length = sizeof(size_t); assert(str_length-1 < UINT16_MAX); // we do p + v + salt[v] in 16-bit elements // TODO: assert((str_length-1)/salt_length * p_increment < UINT16_MAX); __m128i saltv_u8; __m256i v, saltv; if(salt_length == 4) { v = _mm256_set1_epi64x(0x0003000200010000); // `v%salt_length` is 0 1 2 3 repeating saltv_u8 = _mm_set1_epi32( salt_byval ); saltv = _mm256_cvtepu8_epi16( saltv_u8 ); // salt[v] repeats with the same pattern: expand it to 16b elements with pmovzx } else { assert(salt_length == 8); v = _mm256_cvtepu8_epi16( _mm_set1_epi64x(0x0706050403020100) ); saltv_u8 = _mm_set1_epi64x( salt_byval ); saltv = _mm256_cvtepu8_epi16( saltv_u8 ); } __m256i v_saltv = _mm256_add_epi16(v, saltv); __m256i p_increment; __m256i p = init_p_and_increment(salt_length, &p_increment, saltv, saltv_u8); for (unsigned i=str_length-1; i>0 ; /*i-=16 */){ // 16 uint16_t j values per iteration. i-- happens inside the scalar shuffle loop. p = _mm256_add_epi16(p, p_increment); // p += salt[v]; with serial dependencies accounted for, prefix-sum style __m256i j_unmodded = _mm256_add_epi16(v_saltv, p); // size_t j = (v + saltv[v] + p) % i; //////// scalar loop over 16 j elements, doing the modulo and swap // alignas(32) uint16_t jbuf[16]; // portable C++11 syntax uint16_t jbuf[16] __atsortingbute__((aligned(32))); // GNU C syntax _mm256_store_si256((__m256i*)jbuf, j_unmodded); const int jcount = sizeof(jbuf)/sizeof(jbuf[0]); for (int elem = 0 ; elem < jcount ; elem++) { if (--i == 0) break; // in fact returns from the whole function. // 32-bit division is significantly faster than 64-bit division unsigned j = jbuf[elem] % (uint32_t)i; // doubtful that vectorizing this with Intel SVML _mm_rem_epu16 would be a win // since there's no hardware support for it. Until AVX512CD, we need each element in a gp reg as an array index anyway. char temp = str[i]; str[i] = str[j]; str[j] = temp; } } } 

Cela comstack pour asm qui a l'air correct, mais je ne l'ai pas exécuté.

Clang fait une boucle intérieure assez sensible. Ceci est avec -fno-unroll-loops pour la lisibilité. Laissez cela de côté pour la performance, bien que ce ne soit pas important ici car la surcharge de la boucle n'est pas le goulot d'étranglement.

  # The loop part of clang3.8.1's output. -march=haswell -fno-unroll-loops (only for human readability. Normally it unrolls by 2). .LBB0_6: # outer loop # in Loop: Header=BB0_3 Depth=1 add esi, 1 .LBB0_3: # first iteration entry point # =>This Loop Header: Depth=1 vpaddw ymm2, ymm2, ymm1 # p += p_increment vpaddw ymm3, ymm0, ymm2 # v+salt[v] + p vmovdqa ymmword ptr [rsp], ymm3 # store jbuf add esi, -1 lea r8, [rdi + rsi] mov ecx, 1 .LBB0_4: # inner loop # Parent Loop BB0_3 Depth=1 # gcc's version fully unrolls the inner loop, leading to code bloat test esi, esi # if(i==0) return je .LBB0_8 movzx eax, word ptr [rsp + 2*rcx - 2] # load jbuf xor edx, edx div esi mov r9b, byte ptr [r8] # swap mov al, byte ptr [rdi + rdx] mov byte ptr [r8], al mov byte ptr [rdi + rdx], r9b add esi, -1 add r8, -1 cmp rcx, 16 # silly clang, not macro-fusing cmp/jl because it wants to use a weird way to increment. lea rcx, [rcx + 1] jl .LBB0_4 # inner loop jmp .LBB0_6 # outer loop 

[disclaimer: considérant une application 32 bits – dans lequel size_t est unsigned int]

L’allocation alignée peut être effectuée à l’aide de fonctions alignées_malloc, que vous pouvez trouver aussi bien pour Windows que pour Linux .

En affectant votre chaîne de cette manière (à une limite de 64 octets), vous pourrez charger des données directement dans les registres _mm256i en utilisant des charges alignées _mm256_load_si256 pour tous les octets.

Si la chaîne n’est pas correctement alignée, vous pouvez charger les premiers octets à l’aide de charges non alignées _mm256_loadu_si256 pour les octets au début.

La première opération modulo que vous effectuez (v% = salt_length) est effectuée avec un opérande constant que vous pouvez initialiser avant la boucle dans un registre avx à l’aide de _mm256_set1_epi32 :

 __m256i mod = _mm256_set2_epi32(salt_length); 

pour le suivant, vous pouvez utiliser _mm256_set_epi32 qui initialise un registre avec toutes les valeurs fournies (méfiez-vous de l’ordre inverse ).

Notez que si vous utilisez AVX2 ou AVX512 (et pas seulement AVX – votre question est un peu déroutante au sujet du jeu d’instructions), vous pouvez également charger des données à l’aide d’ instructions de regroupement qui chargent les données d’un vecteur à des index spécifiés dans une seconde. argument.

Les plus petits nombres que vous pouvez utiliser avec AVX512 dans les opérations restantes sont en 8 bits, la valeur insortingnsèque étant: _mm512_rem_epu8

Cependant, pour initialiser son entrée avec des valeurs de 0 à N, vous devez utiliser _mm512_set_epi32 et lui transmettre des entiers de 8 bits _mm512_set_epi32 entiers de 32 bits, car il ne semble pas exister insortingnsèquement prenant 64 entiers de 8 bits. Le code ressemblerait à ceci:

 const __m512i sseConst = _mm512_set_epi32( (63<<24) | (62<<16) | (61<<8) | 60, (59<<24) | (58<<16) | (57<<8) | 56, ... etc ...); 

Envisagez d'écrire un générateur de code pour cela, si vous n'aimez pas ce type de frappe ou si vous avez peur des fautes de frappe.

L'utilisation de type __m512i devrait s'occuper automatiquement de l'alignement, si vous ne l' __m512i pas avec malloc() . Sinon, recherchez "aligné malloc" pour votre compilateur. L'alignement dont vous avez besoin est de 64 octets (égal à 512 bits).

Lorsque vous avez besoin des N prochains nombres entiers du vecteur, vous pouvez faire:

 const __m512i inc = _mm512_set1_epi32((N<<24) | (N<<16) | (N<<8) | N); 

Et ensuite, vous pouvez append inc et sseConst ( _mm512_add_epi32 insortingnsèque).