SIMD：累积相邻对

Question

我正在学习如何使用 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 级别，并且该命令会导致下一个级别的精度损失。）

感谢大家。:)

Answer 1

窄类型T=uint8_t或 的uint16_t加宽情况可能最好使用乘数为 的SSSE3pmaddubsw或 SSE2来实现。（内在指南）这些指令是单微操作，并且比洗牌更有效地精确执行您需要的水平加宽添加。pmaddwd1

如果要对 8 个或更多 u8 元素进行水平求和，请psadbw对归零向量使用（绝对差之和）。 _mm_sad_epu8获取每个 64 位元素底部的总和。（如果对整个向量求和，或者对整个 u8 数组求和而不溢出，这是一个很好的第一步。）

如果可以在不损失精度的情况下执行此操作，请先在行之间进行垂直加法，然后再加宽水平加法。（例如，10、12 或 14 位像素分量[u]int16_t不能溢出）。在大多数 CPU 上，负载和垂直相加的每时钟吞吐量为 2（或更多），而pmadd*在 Skylake 及更高版本上，每时钟吞吐量为 1。 这意味着您只需要 1x add + 1x pmadd 与 2x pmadd + 1x add 相比，因此即使在 Skylake 上也是一个重大胜利。 （对于第二种方式，如果你有 AVX，两个加载都可以折叠到 pmadd 的内存操作数中。对于 pmadd 之前的 add 方式，你首先需要一个纯加载，然后将第二个加载折叠到 add 中，所以你可能不会保存前端微指令，除非您使用索引寻址模式并且它们是非层压的。）

理想情况下，您不需要+=进入累加器数组，而只需并行读取 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带有or的 T=U 情况uint32_t_mm_hadd_epi16是 SSSE3或的一个用例_mm_hadd_epi32。它的成本与 2 个随机播放 + 一个垂直添加相同，但无论如何您都需要将 2 个输入打包为 1 个输入。

如果您想解决 Haswell 及更高版本上的 shuffle-port 瓶颈，您可以考虑在输入上使用 qword 移位，然后使用shufps（_mm_shuffle_ps+一些转换）将结果混洗在一起。这可能是 Skylake 的胜利（每个时钟移位吞吐量为 2 个），尽管它总共花费了 5 个 uops，而不是 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);
}

Ice Lake 以及后来的 Intel 和 AMD 对于某些 shuffle 的吞吐量优于 1/clock 吞吐量，但并非全部。例如，vphadddIce Lake 或更高版本上的 p1 或 p5 为 2 uop，对于 p015 中的任何一个加上 1 uop。但是vhaddps2 shuffle uops 只能在端口 5 上运行，即使您可以vhaddps使用 2x 进行精确模拟vshufps（对于_MM_SHUFFLE(2,0,2,0)偶数和_MM_SHUFFLE(3,1,3,1)奇数，如另一个答案所示），并且vshufps在 p1/p5 上运行（vunpcklps由于某种原因不像）。AMD CPU 的运行vhaddps效率也低于您的预期。因此，如果您有大量 FP 数据的成对求和，您可能会考虑手动使用 clang 不会“优化”到vhaddps.

但对于vphadddIce Lake / Alder Lake 和 Zen 3 / Zen 4 上的整数来说，情况似乎没问题。至少有 YMM 向量；vphaddd xmm, xmm, xmm根据https://uops.info/ 的测量，Zen 3/4 至少需要额外的 uop 。

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

hadd32_avx2可以每 3 个时钟运行 1 个（在端口 5 上出现洗牌瓶颈）。用于实现循环的加载+存储+额外的add_epi32微指令可以轻松地在其阴影下运行。

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