SIMD：累积相邻对答案

【问题标题】：SIMD: Accumulate Adjacent PairsSIMD：累积相邻对
【发布时间】：2019-07-30 04:11:06
【问题描述】：

我正在学习如何使用 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 内在函数来解决这个问题？我可以控制 accum 的对齐方式，但不能控制通道。

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

谢谢大家。 :)

【问题讨论】：

看起来像 SSSE3 _mm_hadd_epi32 或 _mm_hadd_epi16 的用例 T 是 int16_t 而不是 int。它的成本与 2 个 shuffle + 一个垂直添加相同，但无论如何您都需要将 2 个输入打包为 1。如果您想解决 Intel CPU 上的 shuffle-port 瓶颈，您可以考虑在输入上使用 qword 移位，然后将结果与shufps 混在一起。
哇，太酷了！我一直认为 SIMD 中不可能进行“水平”操作。我明天试试这个。对于它的价值，此操作的主要用例是 uint8_t -> uint16_t
我没有意识到你正在扩大，这完全改变了事情。（此外，您在 Godbolt 上显示 short -> int；您的目标是什么 SSE/AVX 版本？您在 Godbolt 上使用了 -march=native，即 Skylake-AVX512 = AVX512BW）。无论如何，当U 与T 的宽度不同时，_mm_hadd_* 没有用处。您可能希望 pmaddwd 或 pmaddubsw 的乘数为 1 以将水平对添加到更宽的结果中。
如果您知道这是您的目标，您应该始终使用 -march=haswell 或类似名称。这设置了重要的调整选项以及指令集。并且不要使用-march=corei7，这有点无意义/令人困惑，因为它基本上是-march=nehalem（第一代核心i7）。
在 Godbolt 上，您可以像普通人一样使用 #include <stddef.h> 和使用 size_t。请注意，gcc 确实为 uint8_t -> uint16_t 自动矢量化了您的代码。不是特别很好，但确实做到了。

标签： c++ sse simd intrinsics avx

【解决方案1】：

窄类型 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 add + 1x pmadd 与 2x pmadd + 1x add 因此即使在 Skylake 上也是一个重大胜利。（对于第二种方式，两个负载都可以折叠到 pmadd 的内存操作数中，如果您有 AVX。对于在 pmadd 之前添加方式，您需要先进行纯加载，然后将第二次加载折叠到添加中，因此您可能不会保存前端 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 是打包乘法加字到双字，命名方案与punpcklwd 等相同）

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

如果您想解决 Haswell 及更高版本上的 shuffle-port 瓶颈，您可以考虑在输入上使用 qword 移位，然后将结果与 shufps（_mm_shuffle_ps + 一些转换）混在一起。这可能是 Skylake 的胜利（每个时钟移位吞吐量 2 个），即使它花费了更多 5 个总微指令而不是 3 个。它最多可以运行每个输出向量 5/3 个周期，而不是每个向量运行 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 结果。所以用轮班来模仿 hadd 可能是一个更大的胜利。

// 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 个（饱和所有矢量 ALU 端口）。额外的add_epi32 与accum[] 相加将减慢到每个 256 位结果向量最多 7/3 个周期，并且您需要展开（或使用展开的编译器）而不仅仅是瓶颈前端。

hadd32_avx2 可以每 3 个时钟运行 1 个（对于随机播放，端口 5 存在瓶颈）。用于实现循环的加载 + 存储 + 额外 add_epi32 uops 可以轻松运行。

（https://agner.org/optimize/，见https://stackoverflow.com/tags/x86/info）

【讨论】：