c - 如何初始化范围从0到N的SIMD vector ？

Question

我正在尝试为以下功能编写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。
如何正确初始化向量？

Answer 1

仅仅问如何初始化向量，基本上就是在寻求教程。 Google提出了一些Intel使用内在函数的指南。我将假装这个问题并非微不足道，回答了有关尝试有效地实现此功能的问题。绝对不是您作为初学者矢量化的那种功能。

请参阅x86标签Wiki，以获得指向文档(尤其是esp)的链接。英特尔的内在函数指南。请参阅sse标签Wiki，以获取very nice intro to SIMD programming with SSE intrinsics, and how to use SIMD effectively的链接以及其他链接。

内容摘要:

展开/向量化免费处理v % salt_length。

如果它不是2的幂或编译时常数，如何将v++; v %= loop_invariant;向量化。包括有关标题问题的答案，有关使用_mm_set_epi8或为此目的初始化向量的其他方式。

如何开始矢量化一个复杂的循环，如下所示:从展开查找串行依赖关系开始。

是完整功能的未经测试的版本，可对除% i和swap之外的所有内容进行矢量化处理。 (即，将您便宜的所有操作矢量化，就像您问的那样)。

(v + salt[v] + p)(以及导致它的所有内容)矢量化为两个vpaddw指令。循环外用于对p进行矢量化的前缀总和设置非常棘手，但我最终也对其进行了矢量化。

函数的绝大多数运行时间将在j元素矢量的标量内部循环中，在div上出现瓶颈(或SVML可以执行的任何操作)和/或超大型字符串的高速缓存未命中。

整个循环无法轻易向量化，因为与伪随机索引的交换会创建不可预测的串行依赖关系。可以使用AVX512聚集+随机+分散，并使用AVX512CD查找冲突位掩码，但这必须是一个单独的问题。我不确定要高效地执行此操作有多困难，或者不确定是否经常要重复多次向量随机播放，而只能在一个不冲突的元素中取得进展。

由于 salt_length = sizeof(size_t)是一个编译时常量，并且比向量长度小2的幂，因此 v++和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:，即@harold在揭示关于salt_length的关键信息之前的建议:

由于我们知道 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

即使对于相当大的 p值，这也会迅速使 i太大而无法重复进行 i的条件减法(直到条件掩码全为false)。

有一些通过已知整数常数进行模运算的技巧(请参阅 http://libdivide.com/)，但是AFAIK为附近的除数(即使具有16的2的幂数)计算乘法模逆并不比完全独立的数容易。因此，我们不能便宜地仅调整 i值的下一个向量的常量。

小数定律也许值得剥下最后一对
向量迭代，带有乘法模逆的预先计算的向量
因此 % i可以使用向量完成。一旦我们接近字符串的末尾，它在L1缓存中可能很热，因此我们完全陷入 div的瓶颈，而不是交换加载/存储。为此，我们可能会使用一个序言来获得 i值，该值是16的倍数，因此当我们接近i = 0时，最后几个向量始终具有相同的 i值对齐方式。否则，我们将有一个常数LUT用于一系列i值，并从中简单地进行未对齐的载荷。这意味着我们不必旋转 salt_v和 p。

可能转换为FP很有用，因为最近的Intel CPU(尤其是Skylake)具有功能非常强大的FP分区硬件，并且具有显着的流水线化(吞吐量:延迟率)。如果我们能够通过正确的舍入选择获得准确的结果，那就太好了。 ( float和 double可以精确地表示任何最大为尾数大小的整数。)

我猜值得尝试使用Intel的 _mm_rem_epu16(使用 i的向量递减的 set1(16)值的向量)。如果他们使用FP获得准确的结果，那就太好了。如果仅将其解压缩为标量并进行整数除法，则会浪费时间将值重新返回到向量中。

无论如何，当然，最简单的解决方案是使用标量循环迭代矢量元素。直到您想出使用AVX512CD进行交换的极好的方法之前，这似乎是合理的，但如果它们在L1缓存中都很热，它的速度可能会比交换慢大约一个数量级。

(直到)函数的部分矢量化版本:

这是 the code on the Godbolt compiler explorer，带有完整的设计注释注释，包括我在计算SIMD前缀和算法时所做的图。我最终想起，我曾在 @ZBoson's floating point SSE Prefix sum answer中将它的较窄版本视为构建基块，但直到我自己对其进行了重新发明后，我才意识到。

// See the godbolt link for full design-notes comments
// comments about what makes nice ASM or not.
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>
#include <assert.h>

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，但是我还没有运行它。

叮当作一个相当明智的内循环。这与 -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