SIMD:累积相邻对

SIMD: Accumulate Adjacent Pairs

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

谢谢大家。 :)

窄类型 T = uint8_tuint16_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_tuint32_t 的 T=U 案例是 SSSE3 _mm_hadd_epi16_mm_hadd_epi32 的一个用例.它的成本与 2 次随机播放 + 垂直添加相同,但无论如何您都需要将 2 个输入打包为 1 个。

如果您想解决 Haswell 及更高版本上的洗牌端口瓶颈,您可以考虑对输入使用 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 结果。所以用 shifts 模拟 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 微指令来实现你的循环可以 运行 很容易。

(https://agner.org/optimize/, and see https://whosebug.com/tags/x86/info)