使用 SIMD 找出两个元素的最大差异

Posted 2023-02-16

【中文标题】使用 SIMD 找出两个元素的最大差异

【英文标题】：Using SIMD to find the biggest difference of two elements

【发布时间】：2018-09-24 19:32:35

【问题描述】：

我编写了一个算法来获取 std::vector 中两个元素之间的最大差异，其中两个值中的较大者必须位于比较低值更高的索引处。

    unsigned short int min = input.front();
    unsigned short res = 0;

    for (size_t i = 1; i < input.size(); ++i)
    
        if (input[i] <= min)
        
            min = input[i];
            continue;
        

        int dif = input[i] - min;
        res = dif > res ? dif : res;
    

    return res != 0 ? res : -1;

是否可以使用 SIMD 优化此算法？我是 SIMD 的新手，到目前为止我还没有成功

【参考方案1】：

您没有指定任何特定的架构，因此我将使用英文描述的算法来保持这个大部分架构中立。但它需要一个 SIMD ISA，它可以有效地根据 SIMD 比较结果进行分支，以检查通常为真的条件，例如 x86，但不是真正的 ARM NEON。

这不适用于 NEON，因为它没有等效的移动掩码，而且 SIMD -> 整数会导致许多 ARM 微架构停止。

---

循环遍历数组时的正常情况是一个元素或整个 SIMD 元素向量，不是新的min，而不是diff 候选强>。我们可以快速浏览这些元素，只有在有新的min 时才会放慢速度以获取详细信息。这就像 SIMD strlen 或 SIMD memcmp，除了不是在第一次搜索命中时停止，我们只是对一个块进行标量然后继续。

---

对于输入数组的每个向量 v[0..7]（假设每个向量有 8 个 int16_t 元素（16 字节），但这是任意的）：

SIMD 比较 vmin > v[0..7]，并检查所有元素是否为真。 （例如 x86 _mm_cmpgt_epi16 / if(_mm_movemask_epi8(cmp) != 0)）如果某处有新的min，我们有一个特殊情况：旧的最小值适用于某些元素，但新的最小值适用于其他元素。并且向量中可能有多个新的最小更新，并且在这些点中的任何一个点都有新的差异候选。

所以用标量代码处理这个向量（更新一个标量 diff，它不需要与向量 diffmax 同步，因为我们不需要位置）。

完成后将最终的min 广播到vmin。或者做一个 SIMD 水平 min，这样以后的 SIMD 迭代的乱序执行就可以开始，而无需等待来自标量的 vmin。如果标量代码是无分支的，则应该可以正常工作，因此标量代码中不会出现导致后续向量工作被抛出的错误预测。

作为替代方案，SIMD 前缀和类型的事物（实际上是前缀最小值）可以生成 vmin，其中每个元素都是该点之前的最小值。 (parallel prefix (cumulative) sum with SSE) .您可以总是这样做以避免任何分支，但如果新的最小候选人很少，那么它很昂贵。不过，它在难以分支的 ARM NEON 上可能是可行的。

如果没有新的最小值，SIMD 压缩最大值 diffmax[0..7] = max(diffmax[0..7], v[0..7]-vmin)。 （如果您使用无符号最大值来处理整个范围，请使用饱和减法，这样您就不会得到较大的无符号差异。）

在循环结束时，对diffmax 向量执行 SIMD 水平最大值。请注意，由于我们不需要 需要最大差异的位置，因此我们不需要在找到新的候选者时更新循环内的所有元素。我们甚至不需要保持标量特例 diffmax 和 SIMD vdiffmax 彼此同步，只需在最后检查以获取标量和 SIMD 最大差异的最大值。

---

SIMD min/max 与水平求和基本相同，只是您使用 packed-max 而不是 packed-add。对于 x86，请参阅Fastest way to do horizontal float vector sum on x86。

或者在带有 SSE4.1 的 x86 上用于 16 位整数元素，phminposuw / _mm_minpos_epu16 可用于最小值或最大值、有符号或无符号，并对输入进行适当调整。 max = -min(-diffmax)。您可以将 diffmax 视为无符号，因为它已知为非负数，但 Horizontal minimum and maximum using SSE 显示了如何翻转符号位以将有符号范围移位为无符号并返回。

---

每次我们找到一个新的min 候选者时，我们可能会得到一个分支错误预测，或者我们发现新的min 候选者过于频繁，以至于效率不高。

如果经常期望新的min 候选者，使用较短的向量可能会很好。或者在发现当前向量中有一个新的-min，然后使用更窄的向量仅在更少的元素上进行标量。在 x86 上，您可以使用 bsf（向前位扫描）来查找哪个元素具有第一个 new-min。这使您的标量代码对向量比较掩码具有数据依赖性，但如果错误预测到它的分支，则比较掩码将准备就绪。否则，如果分支预测能够以某种方式找到向量需要标量回退的模式，则预测+推测执行将破坏该数据依赖性。

---

未完成/损坏（由我）示例改编自 @harold 已删除的完全无分支版本的答案，该版本为 x86 SSE2 动态构建了一个最小到该元素的向量。

（@harold 用 suffix-max 而不是 min 写的，我想这就是他删除它的原因。我将它从 max 部分转换为 min。）

x86 的无分支内在函数版本可能看起来像这样某种东西。但是，除非您预期某种斜率或趋势会使新的 min 值频繁出现，否则分支可能会更好。

// BROKEN, see FIXME comments.
// converted from @harold's suffix-max version

int broken_unfinished_maxDiffSSE(const std::vector<uint16_t> &input) 
    const uint16_t *ptr = input.data();

    // construct suffix-min
    // find max-diff at the same time
    __m128i min = _mm_set_epi32(-1);
    __m128i maxdiff = _mm_setzero_si128();

    size_t i = input.size();
    for (; i >= 8; i -= 8) 
        __m128i data = _mm_loadu_si128((const __m128i*)(ptr + i - 8));

   // FIXME: need to shift in 0xFFFF, not 0, for min.
   // or keep the old data, maybe with _mm_alignr_epi8
        __m128i d = data;
        // link with suffix
        d = _mm_min_epu16(d, _mm_slli_si128(max, 14));
        // do suffix-min within block.
        d = _mm_min_epu16(d, _mm_srli_si128(d, 2));
        d = _mm_min_epu16(d, _mm_shuffle_epi32(d, 0xFA));
        d = _mm_min_epu16(d, _mm_shuffle_epi32(d, 0xEE));
        max = d;

        // update max-diff
        __m128i diff = _mm_subs_epu16(data, min);  // with saturation to 0
        maxdiff = _mm_max_epu16(maxdiff, diff);
    

    // horizontal max
    maxdiff = _mm_max_epu16(maxdiff, _mm_srli_si128(maxdiff, 2));
    maxdiff = _mm_max_epu16(maxdiff, _mm_shuffle_epi32(maxdiff, 0xFA));
    maxdiff = _mm_max_epu16(maxdiff, _mm_shuffle_epi32(maxdiff, 0xEE));
    int res = _mm_cvtsi128_si32(maxdiff) & 0xFFFF;

    unsigned scalarmin = _mm_extract_epi16(min, 7);  // last element of last vector
    for (; i != 0; i--) 
        scalarmin = std::min(scalarmin, ptr[i - 1]);
        res = std::max(res, ptr[i - 1] - scalarmin);
    

    return res != 0 ? res : -1;

如果我们处理最后一个完整向量 min 之间的重叠，我们可以用最终未对齐向量替换标量清理。

【讨论】：

那样我做了那个“带后缀的链接”步骤，通过切换到min 被破坏无论如何都是一个不好的方法，真的应该广播然后把它放在 after 块内的东西，使循环携带依赖更快

@harold：哦，好点！这甚至不是 shuffle 的瓶颈，而是 shuffle->max->... dep 链上的瓶颈。我想包含某种代码作为内在函数外观的示例，以便 OP 谷歌更多。因此，我将把那个不太好的代码块留在那里，并带有警告 cmets。如果您有时间改进/修复它，请随时编辑我的答案，或者最好取消删除您自己的答案。