SIMD 搜索最后一个峰值后的低谷

Question

我需要找到比最后一个滚动最大峰值低 X 或更多 % 的值的索引。

峰值是一个数组 ( highs) 中元素的滚动最大值，而值位于另一个数组 ( lows) 中。数组具有相等的长度，并且值保证 <= peaks 数组中的对应项，没有 0、NAN 或无穷大元素。since保证小于till。

迭代实现很简单：

inline
size_t trail_max_intern(double *highs,
        double *lows,
        double max,
        double trail,
        size_t since,
        size_t till)
{
    for (; since < till; ++since) {
        if (max < highs[since]) {
            max = highs[since];
        }

        if (lows[since] / max <= trail) {
            break;
        }
    }

    return since;
}

size_t trail_max_iter(double *highs, double *lows, double trail, size_t since, size_t till)
{
    trail = 1 - trail;
    return trail_max_intern(highs, lows, highs[since], trail, since, till);
}

这基本上是一个线性搜索，条件稍有修改。由于任务上下文（任意的since、till 和trail 值），不能使用其他结构或算法。

为了获得一些加速，我想使用 AVX2 矢量化扩展，看看会发生什么。我的结果是这样的：

size_t trail_max_vec(double *highs, double *lows, double trail, size_t since, size_t till)
{
    double max = highs[since];
    trail = 1 - trail;

    if (till - since > 4) {
        __m256d maxv = _mm256_set1_pd(max);
        __m256d trailv = _mm256_set1_pd(trail);

        for (size_t last = till & ~3; since < last; since += 4) {
            __m256d chunk = _mm256_loadu_pd(highs + since); // load peak block

            // peak rolling maximum computation
            maxv = _mm256_max_pd(maxv, chunk);
            // propagating the maximum value to places 2, 3, and 4
            maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b10010000));
            maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b01000000));
            maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b00000000));

            // divide lows by rolling maximum
            __m256d res = _mm256_div_pd(_mm256_loadu_pd(lows + since), maxv);
            // and if it is lower than the given fraction return its index
            int found = _mm256_movemask_pd(_mm256_cmp_pd(res, trailv, _CMP_LE_OQ));
            if (found) {
                return since + __builtin_ctz(found);
            }

            maxv = _mm256_set1_pd(maxv[3]);
        }

        max = maxv[3];
    }

    // make sure trailing elements are seen
    return trail_max_intern(highs, lows, max, trail, since, till);
}

它产生正确的结果，但比迭代版本慢约 2 倍。我肯定做错了什么，但我不知道是什么。

所以我的问题是，我的方法有什么问题，我该如何解决？

PS 完整源代码可在https://godbolt.org/z/e5YrTo 获得

Answer 1

这会创建一个非常长的依赖链：

        // peak rolling maximum computation
        maxv = _mm256_max_pd(maxv, chunk);
        // propagating the maximum value to places 2, 3, and 4
        maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b10010000));
        maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b01000000));
        maxv = _mm256_max_pd(maxv, _mm256_permute4x64_pd(maxv, 0b00000000));

        // ...

        maxv = _mm256_set1_pd(maxv[3]);

即，您有 8 条指令，全部取决于上一条指令的结果。假设每个都有 3 个周期的延迟，仅此一项就需要 4 个元素的 24 个周期（其余操作可能会在该时间内发生，特别是如果您进行比较lows[since] <= trail * max- 假设max > 0）。

为了减少依赖链，您应该首先计算“本地”滚动最大值，chunk然后计算最大值maxv：

        chunk = _mm256_max_pd(chunk, _mm256_movedup_pd(chunk));                  // [c0, c0c1, c2, c2c3]
        chunk = _mm256_max_pd(chunk, _mm256_permute4x64_pd(chunk, 0b01010100));  // [c0, c0c1, c0c1c2, c0c1c2c3]

        __m256d max_local = _mm256_max_pd(maxv, chunk);

要计算下一个，maxv您可以广播max_local[3]（总延迟约为 6 个周期），或者首先广播chunk[3]并计算其中的最大值maxv（这将只maxpd在依赖链中留下一个，并且您显然会受到吞吐量的限制）。在 Godbolt 上对此进行基准测试会产生很大的噪音，以决定哪个更适合您的情况。

此外，您可以考虑处理更大的块（即加载两个连续的__m256ds，为这些计算区域设置滚动最大值等）