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

Question

我有8个AVX向量，每个向量包含8个浮点(总共64个浮点(，我想将每个向量中的元素求和(基本上执行八个水平求和(。

现在，我使用以下代码：

__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执行通道间混洗。

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

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

Answer 2

Witek902的解决方案应该运行良好，但它可能如果HorizontalSums经常被周围的代码调用，则会遭受高端口5压力。

在Intel Haswell或更新版本上，vhaddps指令解码为3个微操作：2个端口5(p5(微操作和p1或p01的一个微运算(参见Agner Fog的指令表(。函数sort_of_alternative_hadd_ps也解码为3个微操作，但其中只有一个(shuffle(必须在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()内部通过CCD_ 5函数进行回答。总共需要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和CCD_ 7被CPU解码为21个微操作，分别具有13个p5微操作和9个p5微观操作。

根据周围的代码和实际的微体系结构，这种减小的端口5压力可以提高性能。