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

Posted 2023-02-19

【中文标题】计算 8 个 AVX 单精度浮点向量的 8 个水平和

【英文标题】：Computing 8 horizontal sums of eight AVX single-precision floating-point vectors

【发布时间】：2018-12-18 20:01:12

【问题描述】：

我有 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 来执行通道间改组。

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

【问题讨论】：

相关：Most efficient way to get a __m256 of horizontal sums of 8 source __m256 vectors

【参考方案1】：

好的，我想我找到了基于（通常很慢）HADD 的更快算法：

__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 个周期。

【讨论】：

是的，transpose+add 是 HADD 真正成功的用例之一。在我看来很好;你肯定需要在某个地方进行一次过道洗牌，所以我认为你不能避免_mm256_permute2f128_ps 或用vinsertf128 替换它。 （vperm2f128 在 Ryzen 上很慢，但在 Intel 上仍然只有 1 uop。如果为 Ryzen 进行调整，您可能只需使用 128 位向量来减少转置工作量，除非在寄存器中只保存一半的数据是问题。或者对于 Ryzen，extract + insert 会比 vperm2f128 更快，但在 Intel 上当然更慢。）

也许未来的一些 AMD uarch 会根据当前时间将vperm2f128 解码为不同的微指令，但在 Ryzen 上它始终是 8 微指令：/ 有时你可以在不为英特尔牺牲任何东西的情况下为 Ryzen 编写代码，但是这不是那些时代之一。

【参考方案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);

可以替换 Witek902 中的前 4 个 _mm256_hadd_ps() 内在函数 answer 通过 sort_of_alternative_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 的压力可能会提高性能。