c – 最快的步幅-3收集指令序列是什么？
作者：互联网
问题：
从内存生成32位元素的stride-3集合的最有效序列是什么？
如果内存安排如下：
MEM = R0 G0 B0 R1 G1 B1 R2 G2 B2 R3 G3 B3 ...
我们想获得三个YMM寄存器,其中：
YMM0 = R0 R1 R2 R3 R4 R5 R6 R7
YMM1 = G0 G1 G2 G3 G4 G5 G6 G7
YMM2 = B0 B1 B2 B3 B4 B5 B6 B7
动机和讨论
标量C代码就像
template <typename T>
T Process(const T* Input) {
T Result = 0;
for (int i=0; i < 4096; ++i) {
T R = Input[3*i];
T G = Input[3*i+1];
T B = Input[3*i+2];
Result += some_parallelizable_algorithm<T>(R, G, B);
}
return Result;
}
假设some_parallelizable_algorithm是用内在函数编写的,并且可以尽可能快地实现：
template <typename T>
__m256i some_parallelizable_algorithm(__m256i R, __m256i G, __m256i B);
所以T = int32_t的向量实现可以是这样的：
template <>
int32_t Process<int32_t>(const int32_t* Input) {
__m256i Step = _mm256_set_epi32(0, 1, 2, 3, 4, 5, 6, 7);
__m256i Result = _mm256_setzero_si256();
for (int i=0; i < 4096; i+=8) {
// R = R0 R1 R2 R3 R4 R5 R6 R7
__m256i R = _mm256_i32gather_epi32 (Input+3*i, Step, 3);
// G = G0 G1 G2 G3 G4 G5 G6 G7
__m256i G = _mm256_i32gather_epi32 (Input+3*i+1, Step, 3);
// B = B0 B1 B2 B3 B4 B5 B6 B7
__m256i B = _mm256_i32gather_epi32 (Input+3*i+2, Step, 3);
Result = _mm256_add_epi32 (Result,
some_parallelizable_algorithm<int32_t>(R, G, B));
}
// Here should be the less interesting part:
// Perform a reduction on Result and return the result
}
首先,这可以完成,因为有32位元素的收集指令,但16位元素或8位元素没有.
其次,更重要的是,出于性能原因,应完全避免上述收集指令.使用连续的宽负载并对加载的值进行混洗以获得R,G和B向量可能更有效.
template <>
int32_t Process<int32_t>(const int32_t* Input) {
__m256i Result = _mm256_setzero_si256();
for (int i=0; i < 4096; i+=3) {
__m256i Ld0 = _mm256_lddqu_si256((__m256i*)Input+3*i));
__m256i Ld1 = _mm256_lddqu_si256((__m256i*)Input+3*i+1));
__m256i Ld2 = _mm256_lddqu_si256((__m256i*)Input+3*i+2));
__m256i R = ???
__m256i G = ???
__m256i B = ???
Result = _mm256_add_epi32 (Result,
some_parallelizable_algorithm<int32_t>(R, G, B));
}
// Here should be the less interesting part:
// Perform a reduction on Result and return the result
}
似乎对于power-2 strides(2,4,…),已知使用UNKPCKL / UNKPCKH的方法,但对于stride-3访问,我找不到任何引用.
我有兴趣为T = int32_t,T = int16_t和T = int8_t解决这个问题,但为了保持专注,我们只讨论第一种情况.
解决方法:
This article from Intel描述了如何准确地完成您想要的3×8案例.
那篇文章涵盖了浮箱.如果你想要int32,你需要转换输出,因为没有整数版本的_mm256_shuffle_ps().
逐字复制他们的解决方案：
float *p; // address of first vector
__m128 *m = (__m128*) p;
__m256 m03;
__m256 m14;
__m256 m25;
m03 = _mm256_castps128_ps256(m[0]); // load lower halves
m14 = _mm256_castps128_ps256(m[1]);
m25 = _mm256_castps128_ps256(m[2]);
m03 = _mm256_insertf128_ps(m03 ,m[3],1); // load upper halves
m14 = _mm256_insertf128_ps(m14 ,m[4],1);
m25 = _mm256_insertf128_ps(m25 ,m[5],1);
__m256 xy = _mm256_shuffle_ps(m14, m25, _MM_SHUFFLE( 2,1,3,2)); // upper x's and y's
__m256 yz = _mm256_shuffle_ps(m03, m14, _MM_SHUFFLE( 1,0,2,1)); // lower y's and z's
__m256 x = _mm256_shuffle_ps(m03, xy , _MM_SHUFFLE( 2,0,3,0));
__m256 y = _mm256_shuffle_ps(yz , xy , _MM_SHUFFLE( 3,1,2,0));
__m256 z = _mm256_shuffle_ps(yz , m25, _MM_SHUFFLE( 3,0,3,1));
所以这是11条指令. (6次加载,5次洗牌)
在一般情况下,可以在O(S * log(W))指令中进行S x W转置.哪里：
> S是步伐
> W是SIMD宽度
假设存在2向量置换和半向量插入加载,则公式变为：
(S x W load-permute) <= S * (lg(W) + 1) instructions
忽略reg-reg移动.对于像3 x 4这样的退化情况,可能会做得更好.
这是带AVX512的3 x 16负载转换：(6个负载,3个shuffle,6个混合)
FORCE_INLINE void transpose_f32_16x3_forward_AVX512(
const float T[48],
__m512& r0, __m512& r1, __m512& r2
){
__m512 a0, a1, a2;
// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
// 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
// 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
a0 = _mm512_castps256_ps512(_mm256_loadu_ps(T + 0));
a1 = _mm512_castps256_ps512(_mm256_loadu_ps(T + 8));
a2 = _mm512_castps256_ps512(_mm256_loadu_ps(T + 16));
a0 = _mm512_insertf32x8(a0, ((const __m256*)T)[3], 1);
a1 = _mm512_insertf32x8(a1, ((const __m256*)T)[4], 1);
a2 = _mm512_insertf32x8(a2, ((const __m256*)T)[5], 1);
// 0 1 2 3 4 5 6 7 24 25 26 27 28 29 30 31
// 8 9 10 11 12 13 14 15 32 33 34 35 36 37 38 39
// 16 17 18 19 20 21 22 23 40 41 42 43 44 45 46 47
r0 = _mm512_mask_blend_ps(0xf0f0, a0, a1);
r1 = _mm512_permutex2var_ps(a0, _mm512_setr_epi32( 4, 5, 6, 7, 16, 17, 18, 19, 12, 13, 14, 15, 24, 25, 26, 27), a2);
r2 = _mm512_mask_blend_ps(0xf0f0, a1, a2);
// 0 1 2 3 12 13 14 15 24 25 26 27 36 37 38 39
// 4 5 6 7 16 17 18 19 28 29 30 31 40 41 42 43
// 8 9 10 11 20 21 22 23 32 33 34 35 44 45 46 47
a0 = _mm512_mask_blend_ps(0xcccc, r0, r1);
a1 = _mm512_shuffle_ps(r0, r2, 78);
a2 = _mm512_mask_blend_ps(0xcccc, r1, r2);
// 0 1 6 7 12 13 18 19 24 25 30 31 36 37 42 43
// 2 3 8 9 14 15 20 21 26 27 32 33 38 39 44 45
// 4 5 10 11 16 17 22 23 28 29 34 35 40 41 46 47
r0 = _mm512_mask_blend_ps(0xaaaa, a0, a1);
r1 = _mm512_permutex2var_ps(a0, _mm512_setr_epi32( 1, 16, 3, 18, 5, 20, 7, 22, 9, 24, 11, 26, 13, 28, 15, 30), a2);
r2 = _mm512_mask_blend_ps(0xaaaa, a1, a2);
// 0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45
// 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46
// 2 5 8 11 14 17 20 23 26 29 32 35 38 41 44 47
}
逆3 x 16转置商店将作为练习留给读者.
由于S = 3有些退化,所以该模式并不是微不足道的.但是如果你能看到模式,你就可以将它推广到任何奇数S以及任何2的幂.
标签：c,x86,vectorization,avx2 来源： https://codeday.me/bug/20190928/1829073.html