AVX2中有序数组的高效稳定总和

Posted 2023-02-16

【中文标题】AVX2中有序数组的高效稳定总和

【英文标题】：Efficient stable sum of a sorted array in AVX2

【发布时间】：2017-09-08 15:24:17

【问题描述】：

考虑一个由double 数字组成的排序（升序）数组。为了数值稳定性，应该对数组求和，就像从头到尾迭代一样，将总和累加到某个变量中。

如何使用 AVX2 有效地对其进行矢量化？

我已经研究过这个方法 Fastest way to do horizontal vector sum with AVX instructions ，但是将它扩展到一个数组似乎很棘手（可能需要一些分治方法），同时通过确保之前对小数字求和来保持浮点精度将它们添加到更大的数字中。

澄清 1：我认为应该可以，例如将前 4 项相加，然后将它们添加到接下来 4 项的总和中，等等。我愿意用一些稳定性来换取性能。但我更喜欢一种不会完全破坏稳定性的方法。

说明 2：内存不应成为瓶颈，因为数组位于 L3 缓存中（但不在 L1/L2 缓存中，因为数组的各个部分是从不同的线程填充的）。我不想诉诸 Kahan 求和，因为我认为真正重要的是操作的数量，而 Kahan 求和会增加大约 4 倍。

【问题讨论】：

升序求和不会杀死所有并行化空间吗？由于 FP 算术不是关联的，并且您明确谈论数字稳定性，我相信您需要从第一项到最后一项的序列和。顺便说一句，不是我的 DV。

@MargaretBloom，我已经更新了问题并对此进行了澄清。

嗯，如果您愿意牺牲一些稳定性，那么简单的并行求和就可以了，不是吗？您将 d[0]+d[4]+d[8]+..、d[1]+d[5]+d[9]+...等相加。

@SergeRogatch：假设你有很多数字（或者没有必要使用 SIMD）。在这种情况下，添加的数字不应太远。我认为您实际上应该检查一下我的推荐结果，看看您损失了多少精度。

@SergeRogatch：我明白了。我以为你的数字比几百多得多。但在这种情况下，您可以将最后 50-60 个数字相加，结果将是相同的。

【参考方案1】：

如果您需要精确和并行性，请使用 Kahan 求和或其他误差补偿技术来重新排序总和（到具有多个累加器的 SIMD 向量元素步幅）。

正如Twofold fast summation - Evgeny Latkin 指出的那样，如果您在内存带宽上遇到瓶颈，则错误补偿总和不会比最大性能总和慢多少，因为 CPU 有大量计算吞吐量在简单并行化中未使用总结内存带宽的瓶颈

另请参阅（kahan summation avx 的谷歌搜索结果）

https://github.com/rreusser/summation-algorithms

https://scicomp.stackexchange.com/questions/10869/which-algorithm-is-more-accurate-for-computing-the-sum-of-a-sorted-array-of-numb

---

Re：您的想法：按顺序对 4 个数字组求和可以让您隐藏 FP-add 延迟和标量添加吞吐量的瓶颈。

在向量中进行水平求和需要大量的改组，因此这是一个潜在的瓶颈。您可能会考虑加载a0 a1 a2 a3，然后随机获取a0+a1 x a2+a3 x，然后是(a0+a1) + (a2+a3)。你有一个锐龙，对吧？最后一步只是将vextractf128 降至 128b。这仍然是总共 3 个 ADD uops 和 3 个 shuffle uops，但指令比标量或 128b 向量少。

---

您的想法与 Pairwise Summation 非常相似

你总是会得到一些舍入误差，但添加相似数量的数字可以将其最小化。

另请参阅Simd matmul program gives different numerical results，了解有关成对求和和简单高效 SIMD 的一些 cmets。

添加 4 个连续数字与垂直添加 4 个 SIMD 向量之间的差异可能可以忽略不计。 SIMD 向量在数组中为您提供小步幅（SIMD 向量宽度）。除非数组增长得非常快，否则它们的大小仍然基本相似。

您无需在最后进行横向求和即可获得大部分收益。您可以维护 1 或 2 个 SIMD 向量累加器，同时在添加到主累加器之前使用更多 SIMD 寄存器来汇总短期运行（可能是 4 或 8 个 SIMD 向量）。

事实上，让您的总拆分方式更多（跨 SIMD 向量元素）意味着它不会增长得那么大。因此，它有助于解决您要避免的问题，并且水平求和到单个标量累加器实际上会使事情变得更糟，尤其是对于严格排序的数组。

通过乱序执行，您不需要太多的 tmp 累加器来完成这项工作，并隐藏累加到主累加器中的 FP-add 延迟。您可以将几组累积到一个新的tmp = _mm_load_ps() 向量中并将其添加到总数中，OoO exec 将与这些执行重叠。所以你的主循环不需要很大的展开因子。

但它不应该太小，你不想在添加到主累加器时遇到瓶颈。您想限制 FP-add 吞吐量。 （或者，如果您关心 Broadwell/Haswell，并且您的内存带宽并不完全成为瓶颈，则可以将一些 FMA 与 1.0 乘数混合以利用该吞吐量。）

例如Skylake SIMD FP add 具有 4 个周期的延迟，0.5 个周期的吞吐量，因此您需要至少执行 7 个作为短 dep 链的一部分的添加，以便将每个添加到单个累加器中。最好更多，和/或最好有2个长期累加器，以更好地吸收资源冲突造成的调度气泡。

有关多个累加器的更多信息，请参阅_mm256_fmadd_ps is slower than _mm256_mul_ps + _mm256_add_ps?。

【讨论】：

请参阅问题中的说明 2。

@SergeRogatch：看看我首先链接的那篇论文。可能有比 Kahan 更快的选项，但在数值上比普通的多累加器 SIMD 更好。 Ryzen 的 FP 吞吐量只有 Intel CPU 的一半左右，所以是的，你很容易受到来自 L3 的数据的 FP 吞吐量的限制。某种成对求和可能是好的。通过精心挑选的洗牌，您可能会得到很好的结果。

【参考方案2】：

到目前为止，这是我的解决方案：

double SumVects(const __m256d* pv, size_t n) 
  if(n == 0) return 0.0;
  __m256d sum = pv[0];
  if(n == 1) 
    sum = _mm256_permute4x64_pd(sum, _MM_SHUFFLE(3, 1, 2, 0));
   else 
    for(size_t i=1; i+1 < n; i++) 
      sum = _mm256_hadd_pd(sum, pv[i]);
      sum = _mm256_permute4x64_pd(sum, _MM_SHUFFLE(3, 1, 2, 0));
    
    sum = _mm256_hadd_pd(sum, pv[n-1]);
  
  const __m128d laneSums = _mm_hadd_pd(_mm256_extractf128_pd(sum, 1),
    _mm256_castpd256_pd128(sum));
  return laneSums.m128d_f64[0] + laneSums.m128d_f64[1];

一些解释：先将相邻的double数组项相加，如a[0]+a[1]、a[2]+a[3]等，然后将相邻项的和相加。

【讨论】：

我用谷歌搜索的一些东西提到了“成对求和”。这可能与您在这里所做的相同。由于您需要在hadd 之后再进行一次随机播放，您能否通过手动随机播放来以不同的顺序设置垂直add？也许制作一个 Ryzen 版本，它可以做更多 128b 向量的东西，但仍然有一点 256b，以获得更多的前端 uop 吞吐量？

【参考方案3】：

您想玩的游戏可能会适得其反。尝试通过从您最喜欢的分布中生成一堆 iid 样本，对它们进行排序，然后将“按升序求和”与“按升序对每个车道求和，然后对车道求和”进行比较。

对于 4 条车道和 16 条数据，按车道求和大约 28% 的时间会产生较小的误差，而按递增顺序求和会产生大约 17% 的时间的较小误差；对于 4 条车道和 256 条数据，在 68% 的时间内对车道求和会产生较小的误差，而按递增顺序求和会在 12% 的情况下产生较小的误差。对 lanewise 求和也优于您在自我回答中给出的算法。为此，我在 [0,1] 上使用了均匀分布。

【讨论】：

这是一个有趣的观察，但你能发布源代码吗？不清楚您是如何计算误差的。

您是否与 Kahan summation 或其他东西进行了比较以发现错误？或者您是否为float 执行此操作并使用double 或扩展精度获得无错误总和？还是别的什么？

@SergeRogatch：将一堆doubles 相加的误差是计算结果与精确结果之间的差异。使用 mpfr 之类的库可以直接计算出准确的结果。