计算8个AVX单精度浮点向量的8个水平求和

Question

我有8个AVX向量，每个向量包含8个浮点数（总共64个浮点数），我想将每个向量中的元素相加在一起（基本上执行8个水平求和）。

目前，我正在使用以下代码：

__m256 HorizontalSums(__m256 v0, __m256 v1, __m256 v2, __m256 v3, __m256 v4, __m256 v5, __m256 v6, __m256 v7)
{
    // transpose
    const __m256 t0 = _mm256_unpacklo_ps(v0, v1);
    const __m256 t1 = _mm256_unpackhi_ps(v0, v1);
    const __m256 t2 = _mm256_unpacklo_ps(v2, v3);
    const __m256 t3 = _mm256_unpackhi_ps(v2, v3);
    const __m256 t4 = _mm256_unpacklo_ps(v4, v5);
    const __m256 t5 = _mm256_unpackhi_ps(v4, v5);
    const __m256 t6 = _mm256_unpacklo_ps(v6, v7);
    const __m256 t7 = _mm256_unpackhi_ps(v6, v7);

    __m256 v = _mm256_shuffle_ps(t0, t2, 0x4E);
    const __m256 tt0 = _mm256_blend_ps(t0, v, 0xCC);
    const __m256 tt1 = _mm256_blend_ps(t2, v, 0x33);
    v = _mm256_shuffle_ps(t1, t3, 0x4E);
    const __m256 tt2 = _mm256_blend_ps(t1, v, 0xCC);
    const __m256 tt3 = _mm256_blend_ps(t3, v, 0x33);
    v = _mm256_shuffle_ps(t4, t6, 0x4E);
    const __m256 tt4 = _mm256_blend_ps(t4, v, 0xCC);
    const __m256 tt5 = _mm256_blend_ps(t6, v, 0x33);
    v = _mm256_shuffle_ps(t5, t7, 0x4E);
    const __m256 tt6 = _mm256_blend_ps(t5, v, 0xCC);
    const __m256 tt7 = _mm256_blend_ps(t7, v, 0x33);

    // compute sums
    __m256 sum0 = _mm256_add_ps(_mm256_add_ps(tt0, tt1), _mm256_add_ps(tt2, tt3));
    __m256 sum1 = _mm256_add_ps(_mm256_add_ps(tt4, tt5), _mm256_add_ps(tt6, tt7));
    v0 = _mm256_blend_ps(sum0, sum1, 0xF0);
    v1 = _mm256_permute2f128_ps(sum0, sum1, 0x21); // final inter-lane shuffling
    return _mm256_add_ps(v0, v1);
}

你可以看到，我只是把向量转置并在最后求和元素。这里我已经使用了两个技巧：在可能的情况下用_mm256_blend_ps替换_mm256_shuffle_ps，以减少英特尔CPU上5号端口的压力，并且在最后使用_mm256_permute2f128_ps + _mm256_blend_ps来执行跨通道混洗。

有没有更好（更快）的计算方法？

- Witek902

---

相关：获取8个源__m256向量的水平总和的最有效方法 - Peter Cordes

Answer 1

好的，我认为我已经找到了一个更快的算法，基于（通常很慢的）HADDs：

__m256 HorizontalSums(__m256 v0, __m256 v1, __m256 v2, __m256 v3, __m256 v4, __m256 v5, __m256 v6, __m256 v7)
{
    const __m256 s01 = _mm256_hadd_ps(v0, v1);
    const __m256 s23 = _mm256_hadd_ps(v2, v3);
    const __m256 s45 = _mm256_hadd_ps(v4, v5);
    const __m256 s67 = _mm256_hadd_ps(v6, v7);
    const __m256 s0123 = _mm256_hadd_ps(s01, s23);
    const __m256 s4556 = _mm256_hadd_ps(s45, s67);

    // inter-lane shuffle
    v0 = _mm256_blend_ps(s0123, s4556, 0xF0);
    v1 = _mm256_permute2f128_ps(s0123, s4556, 0x21);

    return _mm256_add_ps(v0, v1);
}

根据IACA的数据，Haswell平台上运行速度可以快8个时钟周期左右。

- Witek902

---

是的，转置加是HADD实际上胜利的用例之一。 对我来说看起来不错； 您肯定需要在某个地方进行一个车道交叉洗牌，因此我认为您无法避免 _mm256_permute2f128_ps 或将其替换为 vinsertf128。 （vperm2f128 在Ryzen上很慢，但在Intel上仍然只有1个uop。 可能如果针对Ryzen进行调整，您将只使用128位向量以减少转置工作量，除非在寄存器中仅保存一半的数据是问题。 或者对于Ryzen，提取+插入比 vperm2f128 更快，但当然在Intel上更慢。） - Peter Cordes

也许未来的 AMD 微架构会根据立即数将 vperm2f128 解码为不同的微操作，但在 Ryzen 上它始终是 8 个微操作 :/ 有时你可以编写适用于 Ryzen 的代码，而不需要为 Intel 牺牲任何东西，但这并不是其中之一。 - Peter Cordes

Answer 2

Witek902的solution应该可以很好地工作，但如果周围的代码经常调用HorizontalSums，则可能会遇到高端口5压力的问题。

在Intel Haswell或更新版本中，vhaddps指令解码为3个微操作：2个端口5（p5）微操作和一个p1或p01的微操作（请参见Agner Fog的指令表）。 sort_of_alternative_hadd_ps函数也解码为3个微操作，但只有其中一个（洗牌）必须在p5上执行。

inline __m256 sort_of_alternative_hadd_ps(__m256 x, __m256 y)
{
    __m256 y_hi_x_lo = _mm256_blend_ps(x, y, 0b11001100);      /* y7 y6 x5 x4 y3 y2 x1 x0 */
    __m256 y_lo_x_hi = _mm256_shuffle_ps(x, y, 0b01001110);    /* y5 y4 x7 x6 y1 y0 x3 x2 */
    return _mm256_add_ps(y_hi_x_lo, y_lo_x_hi);
}

可以通过sort_of_alternative_hadd_ps函数替换Witek902的answer中的前4个_mm256_hadd_ps()内置函数。总共需要8个额外的指令来计算水平求和：

__m256 HorizontalSums_less_p5_pressure(__m256 v0, __m256 v1, __m256 v2, __m256 v3, __m256 v4, __m256 v5, __m256 v6, __m256 v7)
{
    __m256 s01 = sort_of_alternative_hadd_ps(v0, v1);
    __m256 s23 = sort_of_alternative_hadd_ps(v2, v3);
    __m256 s45 = sort_of_alternative_hadd_ps(v4, v5);
    __m256 s67 = sort_of_alternative_hadd_ps(v6, v7);
    __m256 s0123 = _mm256_hadd_ps(s01, s23);
    __m256 s4556 = _mm256_hadd_ps(s45, s67);

    v0 = _mm256_blend_ps(s0123, s4556, 0xF0);
    v1 = _mm256_permute2f128_ps(s0123, s4556, 0x21);
    return _mm256_add_ps(v0, v1);
}

这将编译为：

HorizontalSums_less_p5_pressure:
        vblendps        ymm8, ymm0, ymm1, 204
        vblendps        ymm10, ymm2, ymm3, 204
        vshufps ymm0, ymm0, ymm1, 78
        vblendps        ymm9, ymm4, ymm5, 204
        vblendps        ymm1, ymm6, ymm7, 204
        vshufps ymm2, ymm2, ymm3, 78
        vshufps ymm4, ymm4, ymm5, 78
        vshufps ymm6, ymm6, ymm7, 78
        vaddps  ymm0, ymm8, ymm0
        vaddps  ymm6, ymm6, ymm1
        vaddps  ymm2, ymm10, ymm2
        vaddps  ymm4, ymm9, ymm4
        vhaddps ymm0, ymm0, ymm2
        vhaddps ymm4, ymm4, ymm6
        vblendps        ymm1, ymm0, ymm4, 240
        vperm2f128      ymm0, ymm0, ymm4, 33
        vaddps  ymm0, ymm1, ymm0
        ret

最终，Witek902的HorizontalSums和HorizontalSums_less_p5_pressure都会被CPU解码为21个微操作，分别包含13个p5微操作和9个p5微操作。

根据周围代码和实际微架构，这种减少端口5压力的方法可能会提高性能。

- wim