如何初始化范围从0到N的SIMD向量？
我有以下函数我正在尝试编写一个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; } } 我正在尝试向量化v %= salt_length; 。
我想初始化一个包含从0到str_length的数字的向量，以便使用SVML的_mm_rem_epu64来计算每个循环迭代的v。
如何正确初始化向量？
询问如何初始化向量基本上是要求教程。 谷歌了解英特尔使用内在函数的一些指南。 我将假装这个问题并不简单，并回答有关尝试有效实现此function的问题 。 这绝对不是你试图作为一个初学者进行矢量化的那种function。
有关docs的链接，请参阅x86标记wiki。 英特尔的内在指南。 请参阅sse标签wiki，以获取有关使用SSE内在函数的SIMD编程的非常好的介绍的链接，以及如何有效地使用SIMD以及其他链接。
内容摘要 ：
- 展开/矢量化免费处理
v % salt_length。 - 如何向量化
v++; v %= loop_invariant;v++; v %= loop_invariant;如果它不是2的幂或编译时常量。 包含有关使用_mm_set_epi8或其他为此目的初始化矢量的方法的标题问题的答案。 - 如何开始像这样开始向量化复杂循环：从展开开始查找串行依赖关系。
-
一个未经测试的完整函数版本，它可以矢量化除
% i和交换之外的所有内容。 （即向所有便宜的操作进行矢量化，就像你问的那样）。-
(v + salt[v] + p)（以及导致它的所有内容）矢量vpaddw两个vpaddw指令 。 用于向量化p的循环外的前缀和设置是棘手的，但我最终也将它矢量化了。 -
绝大多数函数的运行时间将在
j元素向量的标量内循环中，div瓶颈（或SVML可以做的任何事情），和/或具有非常大的字符串的高速缓存未命中。
-
整个循环不能轻易地向量化，因为具有伪随机索引的交换会产生不可预测的串行依赖性。 使用AVX512采集+ shuffle + scatter，使用AVX512CD查找冲突位掩码，可能是可能的，但这必须是一个单独的问题。 我不确定有效地做这件事有多难，或者如果你经常多次重复一次矢量改组，只能在一个非冲突的元素中取得进展。
由于salt_length = sizeof(size_t)是编译时常量，并且2的幂小于向量长度，因此v++和v%=salt_length根本不需要循环内的任何代码 ，并且作为一个v%=salt_length免费发生 -有效地展开循环以并行执行多个v值的效果。
（使用依赖于平台的salt大小意味着32位构建将无法处理使用64位salt创建的数据。即使x32 ABI具有32位size_t ，因此更改为uint64_t似乎是有意义的，除非你永远不需要在机器之间共享盐渍哈希。）
在标量循环中， v遵循重复模式0 1 2 3 0 1 2 3 …（或对于64位为0..7）。 在向量代码中，我们可以同时使用32字节向量中的4B元素执行8个v值，或者使用2B元素同时执行16次迭代。
因此v成为循环不变的常数向量。 有趣的是， salt[v]也是如此，因此我们永远不需要在循环内进行任何盐表查找。 实际上，可以为标量和向量预先计算v+salt[v] 。
标量版本应预先计算v+salt[v]并展开4或8，删除LUT查找，以便所有内存/缓存吞吐量可用于实际交换。 编译器可能不会为您执行此操作，因此您可能需要手动展开并编写额外的代码来处理最后奇数个字符串字节。 （如果不展开，你仍然可以预先计算v+salt[v]的查找表，其类型足够宽，不会缠绕）。
即使只是确保在编译时已知salt_length也可以提供更好的代码。 v %= compile_time_constant比div insn便宜，而且当它的功率为2时非常便宜。（它只是变成v &= 7 ）。 如果标量版本可以内联，或者您使用salt_length = sizeof(size_t)而不是将其作为函数arg传递，编译器可能会为您执行此操作。
如果您还不知道salt_length：即在您发布有关salt_length的重要信息之前@harold建议的内容：
由于我们知道v < salt_length开头，我们最多只需要一个v -= salt_length将它包装回正确的范围并保持不变。 这被称为“强度降低”操作，因为减法是比除法更弱（和更便宜）的操作。
// 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;
为了向量化：让我们假装我们所知道的只是salt_length <= 16 ，所以我们可以使用32 uint8_t值的向量。 （我们可以使用pshufb来矢量化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. }
当然，如果我们知道salt_length是2的幂，我们应该只掩盖每个元素中的相关低位：
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;
注意到我们从错误的元素大小开始是当我们意识到一次只渲染一行是一个坏主意时，因为现在我们必须改变我们已编写的所有代码以使用更宽的元素，有时需要不同的策略做同样的事。
当然，你需要从一个粗略的轮廓开始，无论是精神上的还是写下来的，你都需要分别完成每一步。 正是在思考这个问题的过程中，您会看到不同部分如何组合在一起。
对于复杂循环， 有用的第一步可能是尝试手动展开标量循环。 这将有助于找到串行依赖关系，以及通过展开简化的事情。
(stuff) % i ：困难的部分
我们需要足够宽的元素来保持i的最大值，因为i 不是 2的幂，而不是常数，因此模运算需要工作。 任何更广泛的浪费都会削减我们的吞吐量。 如果我们可以对循环的其余部分进行矢量化，那么对于不同的str_length范围，甚至可能值得使用不同版本的函数。 （或者可能使用64b元素循环直到i <= UINT32_MAX，然后循环直到i <= UINT16_MAX等）。 如果您知道不需要处理> 4GiB的字符串，则可以通过仅使用32位数学来加速常见情况。 （即使高位全部为零，64位除法也比32位除法慢）。
实际上，我认为我们需要与最大p一样宽的元素 ，因为它会永远累积（直到它在原始标量代码中包含2 ^ 64）。 与恒定模数不同，我们不能只使用p%=i来控制它，即使模数是分布式的。 (123 % 33) % (33-16) != 123 % (33-16) 。 即使对齐到16也无济于事：12345％32！= 12345％48％32
对于i重复条件减法，这将很快使p太大（直到条件掩码全部为假），即使对于相当大的i值也是如此。
有已知整数常数模数的技巧（参见http://libdivide.com/ ），但是AFAIK计算出附近除数的乘法模数逆（即使是像16这样的二次幂的步长）也不是比完全独立的号码更容易。 所以我们不能廉价地调整i值的下一个向量的常数。
小数定律可能使得利用乘法模乘逆的预先计算向量剥离最后几个向量迭代是值得的，因此% i可以用向量完成。 一旦我们接近字符串的末尾，它可能在L1缓存中很热，所以我们在div上完全是瓶颈，而不是交换加载/存储。 为此，我们可能使用序言来达到16的倍数的i值，因此当我们接近i = 0时的最后几个向量总是具有i值的相同对齐。 或者我们对一系列i值有一个常量的LUT，并且只是对它进行未对齐的加载。 这意味着我们不必旋转salt_v和p 。
可能转换为FP将是有用的，因为最近的英特尔CPU（特别是Skylake）具有非常强大的FP分区硬件，具有显着的流水线（吞吐量：延迟比）。 如果我们能够通过正确的舍入选择获得精确的结果，那就太棒了。 （ float和double可以精确地表示任何大约到其尾数大小的整数。）
我想值得尝试英特尔的_mm_rem_epu16 （带有i值的向量，你用set1(16)的向量递减）。 如果他们使用FP来获得准确的结果，那就太好了。 如果它只是解压缩到标量并进行整数除法，那么将浪费时间将值返回到向量中。
无论如何，当然最简单的解决方案是使用标量循环迭代向量元素。 在您使用AVX512CD进行交换之前，你想出了一些非常奇特的东西，这看起来很合理，但如果它们在L1缓存中都很热，它可能比交换速度快一个数量级。
（untest）函数的部分矢量化版本：
这是Godbolt编译器资源管理器上的代码 ，带有完整的设计注释注释，包括我在计算出SIMD前缀和算法时所做的图表。 我最终记得我在@ ZBoson的浮点SSE前缀总和答案中看到了一个较窄版本的这个构建块，但直到我自己重新发明它之后才这样。
// 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 *restrict 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] __attribute__((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; } } }
这编译为看起来正确的asm，但我还没有运行它。
Clang是一个相当明智的内循环。 为了便于阅读，这是使用-fno-unroll-loops 。 把它留给性能，虽然这里没关系，因为循环开销不是瓶颈。
# 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
[免责声明：考虑32位应用程序 – 其中size_t是unsigned int]
可以使用aligned_malloc函数完成对齐分配，您可以在windows和linux中找到它们。
以这种方式分配字符串（到64字节边界）将允许您通过对所有字节使用对齐的加载_mm256_load_si256将数据直接加载到_mm256i寄存器中。
如果字符串未正确对齐，则可以使用未对齐的加载_mm256_loadu_si256来加载第一个字节，以用于开头的字节。
您执行的第一个模运算（v％= salt_length）是使用常量操作数完成的，您可以使用_mm256_set1_epi32在avx寄存器中的循环之前初始化该操作数：
__m256i mod = _mm256_set2_epi32(salt_length);
对于下一个，您可以使用_mm256_set_epi32初始化一个包含所有值的寄存器（注意逆序 ）。
请注意，如果您正在使用AVX2或AVX512（而不仅仅是AVX – 您的问题对指令集有点混淆），您还可以使用收集指令加载数据，这些收集指令在一秒内指定的索引处加载数据论点。
在余数运算中，AVX512可以使用的最小数字是8位，内在的是： _mm512_rem_epu8
但是，要使用从0到N的值初始化其输入，您需要使用_mm512_set_epi32并将其传递给包含在32位整数中的8位整数，因为似乎没有内部采用每个8位的64位整数。 代码看起来像：
const __m512i sseConst = _mm512_set_epi32( (63<<24) | (62<<16) | (61<<8) | 60, (59<<24) | (58<<16) | (57<<8) | 56, ... etc ...);
考虑为此编写代码生成器，如果您不喜欢这样的打字或害怕拼写错误。
如果不使用malloc()分配__m512i类型，则应自动处理对齐。 否则，为您的编译器搜索“aligned malloc”。 您需要的对齐方式是64个字节（等于512位）。
当你需要向量中的下一个N整数时，你可以这样做：
const __m512i inc = _mm512_set1_epi32((N<<24) | (N<<16) | (N<<8) | N);
然后你可以添加inc和sseConst （内在_mm512_add_epi32 ）。