使用 SIMD 根据另一个向量位值计算值的乘积

Posted 2023-02-16

【中文标题】使用 SIMD 根据另一个向量位值计算值的乘积

【英文标题】：Computing a product of values based on another vector bit values using SIMD

【发布时间】：2019-11-07 17:17:32

【问题描述】：

我有两个向量。大小为 N 的双精度向量 a 和大小为 ceil(N/8) 的无符号字符向量 b。目标是计算a 的一些值的乘积。 b 将逐位读取，其中每个位指示是否要在产品中考虑从a 给出的double。

  // Let's create some data      
  unsigned nbBits  = 1e7;
  unsigned nbBytes = nbBits / 8;
  unsigned char nbBitsInLastByte = nbBits % 8;
  assert(nbBits == nbBytes * 8 + nbBitsInLastByte);
  std::vector<double> a(nbBits, 0.999999);   // In practice a values will vary. It is just an easy to build example I am showing here
  std::vector<unsigned char> b(nbBytes, false); // I am not using `vector<bool>` nor `bitset`. I've got my reasons!
  assert(a.size() == b.size() * 8);

  // Set a few bits to true
  for (unsigned byte = 0 ; byte < (nbBytes-1) ; byte+=2)
  
    b[byte] |= 1 << 2; // set second (zero-based counting) bit to 'true'
    b[byte] |= 1 << 7; // set last bit to 'true'
                //  ^ This is the bit index
  

如上所述，我的目标是每当b 为真时计算a 中的值的乘积。它可以通过

来实现

  // Initialize the variable we want to compute
  double product = 1.0;

  // Product for the first nbByts-1 bytes
  for (unsigned byte = 0 ; byte < (nbBytes-1) ; ++byte)
  
    for (unsigned bit = 0 ; bit < 8 ; ++bit) // inner loop could be manually unrolled
    
      if((b[byte] >> bit) & 1) // gets the bit value
        product *= a[byte*8+bit];
    
  

  // Product for the last byte
  for (unsigned bit = 0 ; bit < nbBitsInLastByte ; ++bit)
  
    if((b[nbBytes-1] >> bit) & 1) // gets the bit value
      product *= a[(nbBytes-1)*8+bit];
  

这个乘积计算是我的代码中比较慢的部分。我想知道显式矢量化（SIMD）过程是否可以在这里有所帮助。我一直在查看“xmmintrin.h”中提供的功能，但我对 SIMD 了解不多，也没有找到有用的东西。你能帮帮我吗？

【问题讨论】：

所以你正在使用mul而不是更常见的add进行蒙版缩减。相关：Fast dot product of a bit vector and a floating point vector，我想我记得最近的一个浮点 * 布尔点积问题（当然也可以归结为浮点向量的蒙面缩减。哦，找到了：Vector matrix multiplication, float vector, binary matrix，比我记得的要少。另见is there an inverse instruction to the movemask instruction in intel avx2?

有了这几个，我会首先考虑在外循环的开头添加一个if(b[byte]==0) continue;（大约 92% 的时间都是这样）。对于!=0 的情况，可能会为每个（半个）字节进行切换，并使用一些手动调整的实现。乘法的延迟可能不是问题，但分支预测是。如果您的位模式是随机的，大约有 50%，请考虑 Peter 链接到的选项之一。

@chtz：乘法延迟可能很重要！现代 x86 在乘法延迟与吞吐量方面有 8 或 10 倍的差异。即使是位图的开销 -> 矢量掩码也不一定会隐藏每次乘法 4 到 5 个周期的延迟。

刚刚意识到像 add 那样简单地屏蔽是行不通的：乘法的标识元素是 1.0 而不是 0.0。所以你可能不得不融入1.0。 （使用 AVX512，将掩蔽合并到目的地仍然可以正常工作）。如果没有 AVX512，您可能希望将 1.0 混合到输入中，而不是通过混合向量的旧值来增加乘法 dep 链的延迟。无论如何，是的，在一组全为零的掩码位上进行分支值得尝试。是的，SSE2 是 x86-64 的基准。 AVX 和 AVX2 很常见，但除了服务器端的东西，你不能假设它

@chtz：前端的分支错误预测惩罚可能是 10 到 15 个周期，在此期间您可以使用 AVX 执行 10 * 2 * 4 双精度乘法。因此，如果您遇到前端瓶颈，您可能会损失大量吞吐量。但是由于后端吞吐量瓶颈，基于此检查的早期退出关键路径可能只会给后端一些时间赶上而不会损失太多吞吐量。 （在 Nehalem 或更高版本上，可以快速恢复分支未命中，而不是等到分支退休）。手波浪形；必须以真实数据为基准

【参考方案1】：

如果您不关心乘法顺序，那很容易。关键是来自 SSE 4.1 集的_mm_blendv_pd 指令。这使您可以完全无分支。

// Load 2 double values from source pointer, and conditionally multiply with the product.
// Returns the new product.
template<int startIdx>
inline __m128d product2( const double* pSource, __m128i mask, __m128d oldProduct )

    // Multiply values unconditionally
    const __m128d source = _mm_loadu_pd( pSource + startIdx );
    const __m128d newProduct = _mm_mul_pd( source, oldProduct );

    // We only calling product2 with 4 different template arguments.
    // There are 16 vector registers in total, enough for all 4 different `maskAndBits` values.
    constexpr int64_t bit1 = 1 << startIdx;
    constexpr int64_t bit2 = 1 << ( startIdx + 1 );
    const __m128i maskAndBits = _mm_setr_epi64x( bit1, bit2 );
    mask = _mm_and_si128( mask, maskAndBits );

    // NAN if the mask is 0 after the above AND i.e. the bit was not set, 0.0 if the bit was set
    const __m128d maskDouble = _mm_castsi128_pd( _mm_cmpeq_epi64( mask, _mm_setzero_si128() ) );

    // This instruction actually does the masking, it's from SSE 4.1
    return _mm_blendv_pd( newProduct, oldProduct, maskDouble );


double conditionalProducts( const double* ptr, const uint8_t* masks, size_t size )

    // Round down the size of your input vector, and multiply last couple values the old way.
    assert( 0 == size % 8 );
    __m128d prod = _mm_set1_pd( 1.0 );
    const double* const end = ptr + size;
    while( ptr < end )
    
        // Broadcast the mask byte into 64-bit integer lanes
        const __m128i mask = _mm_set1_epi64x( *masks );
        // Compute the conditional products of 8 values
        prod = product2<0>( ptr, mask, prod );
        prod = product2<2>( ptr, mask, prod );
        prod = product2<4>( ptr, mask, prod );
        prod = product2<6>( ptr, mask, prod );
        // Advance the pointers
        ptr += 8;
        masks++;
    
    // Multiply two lanes together
    prod = _mm_mul_sd( prod, _mm_shuffle_pd( prod, prod, 0b11 ) );
    return _mm_cvtsd_f64( prod );

【讨论】：

由于bendv_pd 仅使用第64 位和第128 位作为条件，因此我提出了不同的解决方案。应该可以将 uint8_t 掩码转换为 uint64_t，将其移动 2 个不同的量，以便相关位最终位于第 64 位和第 128 位，_mm_setr_epi64x

@Sopel 不确定我是否关注。掩码中有 8 个相关位。您是否建议将值从 GP 寄存器复制到向量寄存器 4 次而不是 1 次？我的代码中的pand 和pcmpeqq 都有1 个周期延迟/每个。

godbolt.org/z/WxG8Cz 不确定是否可以在没有 AVX2 的情况下有效完成，因为没有 sll 每个通道的计数。 const __m128i mulmask = _mm_sllv_epi64(mask, _mm_set_epi64x((63 - startIdx), (63 - ( startIdx + 1 )))); 将相关条件位获取到 64 位通道的最高位 - 由 blendv 进行检查

@Sopel 我不确定您的代码是否更快。 AVX2 指令vpsllvq 在除 Skylake 之外的大多数 CPU 上都不是特别快。它在 Haswell 和 Broadwell 上有 2 个周期延迟，在 Ryzen 上有 3 个周期延迟。吞吐量也不是很大，它解码为 2-3 µops。我的代码中的pand 和pcmpeqq 都解码为单个 µop/each，它们有 1 个周期延迟，并且吞吐量很好，它们每个周期可以运行 2-3 次。位掩码的内容在对 product2 的 4 次调用中没有数据依赖关系，流水线有望完成这项工作。

是的，你可能是对的。我没有对它进行基准测试，只是一个想法。