我有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来执行通道间混洗。
有没有更好(更快)的方法来计算它?
好的,我想我找到了基于(通常很慢)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个周期。
Witek902的solution应该运行良好,但如果周围的代码经常调用HorizontalSums
,它可能会受到高端口压力的影响。
在Intel Haswell或更新版本中,vhaddps
指令解码为3个微操作:2个端口5(p5)微操作和1个微操作用于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);
}
有可能用_mm256_hadd_ps()
函数替换Witek902的answer中的前4个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压力可以改善性能。