c++ - SIMD:累积相邻对

我正在学习如何使用SIMD内部函数和自动向量化。幸运的是，我有一个正在从事的有用项目，对于SIMD来说，它似乎极为适合，但对于像我这样的新手来说仍然很棘手。

我正在为图像编写一个过滤器，以计算2x2像素的平均值。我正在通过将两个像素的总和累加为一个像素来进行计算。

template <typename T, typename U>
inline void accumulate_2x2_x_pass(
  T* channel, U* accum,
  const size_t sx, const size_t sy,
  const size_t osx, const size_t osy,
  const size_t yoff, const size_t oyoff
) {

  const bool odd_x = (sx & 0x01);

  size_t i_idx, o_idx;

  // Should be vectorizable somehow...
  for (size_t x = 0, ox = 0; x < sx - (size_t)odd_x; x += 2, ox++) {
    i_idx = x + yoff;
    o_idx = ox + oyoff;
    accum[o_idx] += channel[i_idx];
    accum[o_idx] += channel[i_idx + 1];
  }

  if (odd_x) {
    // << 1 bc we need to multiply by two on the edge
    // to avoid darkening during render
    accum[(osx - 1) + oyoff] += (U)(channel[(sx - 1) + yoff]) * 2;
  }
}

但是，Godbolt显示我的循环无法自动矢量化。 (https://godbolt.org/z/qZxvof)如何构造SIMD内部函数来解决此问题？我可以控制累积的对齐方式，但不能控制通道。

(我知道这里有一个平均内在函数，但这在这里不合适，因为我需要生成多个mip级别，并且该命令会导致下一级别的精度损失。)

感谢大家。 :)

最佳答案

窄类型T = uint8_t或uint16_t 的扩大情况可能最好与SSSE3 pmaddubsw或SSE2 pmaddwd乘以1的乘数一起实现。 (Intrinsics guide)这些指令是单码执行的，并且精确地进行水平加宽会比混洗更有效。

如果可以这样做，而又不损失精度，请先在行之间进行垂直添加，然后再扩大水平添加。 (例如[u]int16_t中的10、12或14位像素分量不会溢出)。在大多数CPU上，负载和垂直添加的时钟速率至少为2，而pmadd*的时钟速率仅为每个时钟1，而Skylake及更高版本的时钟速率仅为2。 ，这意味着您只需要添加1x + 1x pmadd与2x pmadd + 1x add，即使在Skylake上，也是一个巨大的胜利。 (对于第二种方式，如果您有AVX，则两个负载都可以折叠到pmadd的内存操作数中。对于在pmadd方式之前进行添加，首先需要一个纯负载，然后将第二个负载折叠成add，所以您可能不需要保存前端uops，除非您使用索引寻址模式并且它们取消分层。)

理想情况下，您不需要将+=放入一个累加器数组中，而是可以并行读取2行并且累加器是只写的，因此您的循环只有2个输入流和1个输出流。

// SSSE3
__m128i hadd_widen8_to_16(__m128i a) {
                      // uint8_t, int8_t  (doesn't matter when multiplier is +1)
    return _mm_maddubs_epi16(a, _mm_set_epi8(1));
}

// SSE2
__m128i hadd_widen16_to_32(__m128i a) {
                   // int16_t, int16_t
    return _mm_madd_epi16(a, _mm_set_epi16(1));
}

这些端口直接连接到256位AVX2，因为输入和输出宽度相同。无需改组固定车道内包装。

是的，它们都是_epi16。英特尔可能与内在名称完全不一致。 asm助记符更一致，更容易记住什么。 (ubsw =无符号字节到带符号的字，除了输入之一是带符号的字节。pmaddwd打包后乘以将双字添加到dword，命名方式与punpcklwd相同，等等。)

带有uint16_t或uint32_t的T = U案例是SSSE3 _mm_hadd_epi16或_mm_hadd_epi32的用例。它的成本与2个混洗+垂直添加相同，但是无论如何您都需要将2个输入打包为1个。

如果要解决Haswell及更高版本上的shuffle-port瓶颈问题，可以考虑在输入上使用qword shift，然后将结果与shufps(_mm_shuffle_ps +一些强制转换)一起改组。这可能是在Skylake上的一个胜利(每个时钟移位吞吐量有2个)，即使它花费了5而不是3的总uops。它可以在每个输出 vector 上最好以5/3周期运行，而不是每个 vector 2个周期。没有前端瓶颈

// UNTESTED

//Only any good with AVX, otherwise the extra movdqa instructions kill this
//Only worth considering for Skylake, not Haswell (1/c shifts) or Sandybridge (2/c shuffle)
__m128i hadd32_emulated(__m128i a, __m128i b) {
    __m128i a_shift = _mm_srli_epi64(a, 32);
    __m128i b_shift = _mm_srli_epi64(b, 32);
    a = _mm_add_epi32(a, a_shift);
    b = _mm_add_epi32(b, b_shift);
    __m128 combined = _mm_shuffle_ps(_mm_castsi128_ps(a), _mm_castsi128_ps(b), _MM_SHUFFLE(2,0,2,0));
    return _mm_castps_si128(combined);
}

对于AVX2版本，您需要横穿通道进行混洗以修正vphadd结果。因此，在轮类中效仿哈德可能是一个更大的胜利。

// 3x shuffle 1x add uops
__m256i hadd32_avx2(__m256i a, __m256i b) {
    __m256i hadd = _mm256_hadd_epi32(a, b);  // 2x in-lane hadd
    return _mm256_permutex_epi64( hadd, _MM_SHUFFLE(3,1,2,0) );
}

// UNTESTED
// 2x shift, 2x add, 1x blend-immediate (any ALU port), 1x shuffle
__m256i hadd32_emulated_avx2(__m256i a, __m256i b)
{
        __m256i a_shift = _mm256_srli_epi64(a, 32);  // useful result in the low half of each qword
        __m256i b_shift = _mm256_slli_epi64(b, 32);  // ... high half of each qword
        a = _mm256_add_epi32(a, a_shift);
        b = _mm256_add_epi32(b, b_shift);
        __m256i blended = _mm256_blend_epi32(a,b, 0b10101010);  // alternating low/high results
        return _mm256_permutexvar_epi32(_mm256_set_epi32(7,5,3,1, 6,4,2,0),  blended);
}

在Haswell和Skylake上，hadd32_emulated_avx2可以每2个时钟1个运行(使所有 vector ALU端口饱和)。合计为add_epi32的额外accum[]会将其减慢到每个256位结果 vector 最多7/3个周期，并且您需要展开(或使用展开的编译器)，而不仅仅是前端的瓶颈。
hadd32_avx2可以以每3个时钟1个的频率运行(端口5出现瓶颈，以进行随机播放)。用于实现循环的load + store +额外的add_epi32 uops可以很容易地在其中运行。

(https://agner.org/optimize/，请参阅https://stackoverflow.com/tags/x86/info)

关于c++ - SIMD:累积相邻对，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/55057933/