使用 SIMD 解决循环数据依赖性 - 在 int8_t sgn 值数组中查找 -1 和 +1 之间的转换
Solve loop data dependency with SIMD - finding transitions between -1 and +1 in an int8_t array of sgn values
我尝试实现性能提升并在 SIMD 方面取得了一些不错的体验。到目前为止，我一直在使用 OMP，并希望使用内在函数进一步提高我的技能。
在以下场景中，由于测试元素 n+1 所需的 last_value 的数据依赖性，我未能改进（甚至矢量化）。
环境是具有 AVX2 的 x64，所以想找到一种方法来矢量化和 SIMDfy 这样的函数。
inline static size_t get_indices_branched(size_t* _vResultIndices, size_t _size, const int8_t* _data) {
size_t index = 0;
int8_t last_value = 0;
for (size_t i = 0; i < _size; ++i) {
if ((_data[i] != 0) && (_data[i] != last_value)) {
// add to _vResultIndices
_vResultIndices[index] = i;
last_value = _data[i];
++index;
}
}
return index;
}
输入是一个有符号 1 字节值的数组。每个元素都是 <0,1,-1> 之一。
输出是输入值（或指针）的索引数组，表示更改为 1 或 -1。
示例in/output
in: { 0,0,1,0,1,1,-1,1, 0,-1,-1,1,0,0,1,1, 1,0,-1,0,0,0,0,0, 0,1,1,1,-1,-1,0,0, ... }
out { 2,6,7,9,11,18,25,28, ... }
我的第一次尝试是尝试各种无分支版本，并通过比较汇编输出来查看自动矢量化或 OMP 是否能够将其转换为 SIMDish 代码。
示例尝试
int8_t* rgLast = (int8_t*)alloca((_size + 1) * sizeof(int8_t));
rgLast[0] = 0;
#pragma omp simd safelen(1)
for (size_t i = 0; i < _size; ++i) {
bool b = (_data[i] != 0) & (_data[i] != rgLast[i]);
_vResultIndices[index] = i;
rgLast[i + 1] = (b * _data[i]) + (!b * rgLast[i]);
index += b;
}
由于没有实验产生 SIMD 输出，我开始试验内在函数，目的是将条件部分转换为掩码。
对于 != 0 部分非常简单：
__m256i* vData = (__m256i*)(_data);
__m256i vHasSignal = _mm256_cmpeq_epi8(vData[i], _mm256_set1_epi8(0)); // elmiminate 0's
测试“最后一次翻转”的条件方面我还没有找到方法。
为了解决以下输出打包方面，我假设 可以工作。
更新 1
深入研究这个话题会发现，分离 1/-1 并去掉 0 是有益的。
幸运的是，在我的例子中，我可以直接从预处理中获取它并使用 _mm256_xor_si256 跳过处理到 <1,0,-1>，例如，将 2 个输入向量分隔为 gt0（全 1）和 lt0 （全部为 -1）。这也允许 4 倍更紧密的数据打包。
我可能想以这样的过程结束
现在的挑战是如何根据 gt0 和 lt0 掩码创建过渡掩码。
更新 2
显然，一种将 1 和 -1 分成 2 个流的方法（参见如何回答）在访问元素以交替扫描时引入了依赖性：
使用
创建@aqrit 转换掩码
transition mask = ((~lt + gt) & lt) | ((~gt + lt) & gt) 是可以的。尽管这增加了相当多的指令，但它似乎是消除数据依赖性的有益权衡。我假设寄存器越大增益越大（可能取决于芯片）。
更新 3
通过向量化 transition mask = ((~lt + gt) & lt) | ((~gt + lt) & gt)
我可以编译这个输出
vmovdqu ymm5,ymmword ptr transition_mask[rax]
vmovdqu ymm4,ymm5
vpandn ymm0,ymm5,ymm6
vpaddb ymm1,ymm0,ymm5
vpand ymm3,ymm1,ymm5
vpandn ymm2,ymm5,ymm6
vpaddb ymm0,ymm2,ymm5
vpand ymm1,ymm0,ymm5
vpor ymm3,ymm1,ymm3
vmovdqu ymmword ptr transition_mask[rax],ymm3
乍一看，与 post 处理（垂直扫描 + 附加到输出）的潜在条件相关陷阱相比，它似乎更有效，尽管处理 2 个流而不是 1 个流似乎是正确且合乎逻辑的.
这缺少每个周期生成初始状态的能力（从 0 到 1 或 -1 的转换）。
不确定是否有一种方法可以增强 transition_mask 生成“位旋转”，或者使用 auto initial _tzcnt_u32(mask0) > _tzcnt_u32(mask1)，就像 Soons 在此处使用的那样： 似乎包含一个分支。
结论
@aqrit 分享的方法是使用改进的 bit-twiddling 解决方案为每个块加载找到转换，结果证明是运行时性能最高的方法。 热内循环只有 9 条 asm 指令长（每 2 个找到的项目与其他方法进行比较）使用 tzcnt 和 blsr 像这样
tzcnt rax,rcx
mov qword ptr [rbx+rdx*8],rax
blsr rcx,rcx
tzcnt rax,rcx
mov qword ptr [rbx+rdx*8+8],rax
blsr rcx,rcx
add rdx,2
cmp rdx,r8
jl main+2580h (...)
对于您的情况，完全矢量化不是最佳选择。这在技术上是可行的，但我认为生成 uint64_t 值数组的开销（我假设你正在为 64 位 CPU 编译）会吃掉所有的利润。
相反，您应该加载 32 字节的块，并立即将它们转换为位掩码。方法如下：
inline void loadBits( const int8_t* rsi, uint32_t& lt, uint32_t& gt )
{
const __m256i vec = _mm256_loadu_si256( ( const __m256i* )rsi );
lt = (uint32_t)_mm256_movemask_epi8( vec );
const __m256i cmp = _mm256_cmpgt_epi8( vec, _mm256_setzero_si256() );
gt = (uint32_t)_mm256_movemask_epi8( cmp );
}
您的其余代码应处理这些位图。要找到第一个 non-zero 元素（您只需在数据的开头这样做），扫描 (lt | gt) 整数中的最低有效设置位。要查找 -1 数字，扫描 lt 整数中的最低有效设置位，查找 +1 数字扫描 gt 整数中的最低有效设置位。找到并处理后，您可以使用按位与清除两个整数的低位部分，或者将它们都向右移动。
CPU 有 BSF instruction which scans for the lowest set bit in an integer, and returns two things at once: a flag telling if the integer was zero, and the index of that set bit. If you’re using VC++ there’s _BitScanForward 内在，否则使用内联 ASM，该指令仅在 VC++ 中公开； GCC 的 __builtin_ctz 不是一回事，它只是 returns 一个值而不是两个值。
但是，在 AMD CPU 上，TZCNT instruction from BMI 1 set is somewhat faster than the old-school BSF（在 Intel 上它们是相等的）。在 AMD 上，TZCNT 可能会稍微快一些，尽管有额外的指令与 0 进行比较。
在 64 位 SIMD 通道之间串行传输状态比在 64 位通用寄存器 (gpr) 之间串行传输状态更昂贵。
实际上，lookup-tables（或 SIMD left-packing）一次只能处理 8 个元素。如果数据平均每 64 个元素保留 6 个元素，那么 left-packing 会浪费大量处理
（特别是如果我们正在收集偏移量而不是进行收集操作）。如果位集密集，则考虑移动到 lookup-tables.
正如@Snoots 所建议的那样，使用 SIMD 创建 64 位位集并使用位扫描内在函数查找所需位集的索引。
分支预测错误：
使用 transition_mask = ((~lt + gt) & lt) | ((~gt + lt) & gt) 或@FalkHüffner transition_mask = (lt ^ (lt - gt)) & (gt ^ (gt – lt)).
其中一项算术运算的状态为 carry-in/carry-out。我会小心使用 _subborrow_u64，因为它是相当不常见的内在函数（并且在旧编译器上存在错误）。
剩下唯一的分支在位扫描操作上循环。必须提取所有设置位.. 但我们可以展开操作和超调以使分支更可预测。超调量需要调整到预期的数据集。
未测试。未检查 Asm。
#include <immintrin.h>
#include <stdint.h>
static inline
uint64_t get_mask (int8_t* src, unsigned char* state) {
__m256i src0 = _mm256_loadu_si256((__m256i*)(void*)src);
__m256i src1 = _mm256_loadu_si256((__m256i*)(void*)&src[32]);
uint64_t lt = (uint32_t)_mm256_movemask_epi8(src0) |
(((uint64_t)(uint32_t)_mm256_movemask_epi8(src1)) << 32);
src0 = _mm256_cmpgt_epi8(src0, _mm256_setzero_si256());
src1 = _mm256_cmpgt_epi8(src1, _mm256_setzero_si256());
uint64_t gt = (uint32_t)_mm256_movemask_epi8(src0) |
(((uint64_t)(uint32_t)_mm256_movemask_epi8(src1)) << 32);
// if borrow then greater-than span extends past the msb
uint64_t m;
unsigned char s = *state;
*state = _subborrow_u64(s, lt, gt, (unsigned long long*)&m); // sbb
return (m ^ lt) & ((gt - (lt + !s)) ^ gt);
}
static inline
size_t bitset_to_index (uint64_t* dst, uint64_t base, uint64_t mask) {
int64_t cnt = _mm_popcnt_u64(mask);
int64_t i = 0;
do { // unroll to taste...
dst[i + 0] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 1] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 2] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 3] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
i += 4;
} while (i < cnt);
return (size_t)cnt;
}
static
uint64_t* get_transition_indices (uint64_t* dst, int8_t* src, size_t len) {
unsigned char state = 0; // in less-than span
uint64_t base = 0; // offset into src array
size_t end = len / 64;
for (size_t i = 0; i < end; i++) {
uint64_t mask = get_mask(src, &state);
src += 64;
dst += bitset_to_index(dst, base, mask);
base += 64;
}
if (len % 64) {
; // todo: tail loop
}
return dst;
}
我尝试实现性能提升并在 SIMD 方面取得了一些不错的体验。到目前为止，我一直在使用 OMP，并希望使用内在函数进一步提高我的技能。
在以下场景中，由于测试元素 n+1 所需的 last_value 的数据依赖性，我未能改进（甚至矢量化）。
环境是具有 AVX2 的 x64，所以想找到一种方法来矢量化和 SIMDfy 这样的函数。
inline static size_t get_indices_branched(size_t* _vResultIndices, size_t _size, const int8_t* _data) {
size_t index = 0;
int8_t last_value = 0;
for (size_t i = 0; i < _size; ++i) {
if ((_data[i] != 0) && (_data[i] != last_value)) {
// add to _vResultIndices
_vResultIndices[index] = i;
last_value = _data[i];
++index;
}
}
return index;
}
输入是一个有符号 1 字节值的数组。每个元素都是 <0,1,-1> 之一。 输出是输入值（或指针）的索引数组，表示更改为 1 或 -1。
示例in/output
in: { 0,0,1,0,1,1,-1,1, 0,-1,-1,1,0,0,1,1, 1,0,-1,0,0,0,0,0, 0,1,1,1,-1,-1,0,0, ... }
out { 2,6,7,9,11,18,25,28, ... }
我的第一次尝试是尝试各种无分支版本，并通过比较汇编输出来查看自动矢量化或 OMP 是否能够将其转换为 SIMDish 代码。
示例尝试
int8_t* rgLast = (int8_t*)alloca((_size + 1) * sizeof(int8_t));
rgLast[0] = 0;
#pragma omp simd safelen(1)
for (size_t i = 0; i < _size; ++i) {
bool b = (_data[i] != 0) & (_data[i] != rgLast[i]);
_vResultIndices[index] = i;
rgLast[i + 1] = (b * _data[i]) + (!b * rgLast[i]);
index += b;
}
由于没有实验产生 SIMD 输出，我开始试验内在函数，目的是将条件部分转换为掩码。
对于 != 0 部分非常简单：
__m256i* vData = (__m256i*)(_data);
__m256i vHasSignal = _mm256_cmpeq_epi8(vData[i], _mm256_set1_epi8(0)); // elmiminate 0's
测试“最后一次翻转”的条件方面我还没有找到方法。
为了解决以下输出打包方面，我假设
更新 1
深入研究这个话题会发现，分离 1/-1 并去掉 0 是有益的。
幸运的是，在我的例子中，我可以直接从预处理中获取它并使用 _mm256_xor_si256 跳过处理到 <1,0,-1>，例如，将 2 个输入向量分隔为 gt0（全 1）和 lt0 （全部为 -1）。这也允许 4 倍更紧密的数据打包。
我可能想以这样的过程结束
更新 2
显然，一种将 1 和 -1 分成 2 个流的方法（参见如何回答）在访问元素以交替扫描时引入了依赖性：
使用
创建@aqrit 转换掩码
transition mask = ((~lt + gt) & lt) | ((~gt + lt) & gt) 是可以的。尽管这增加了相当多的指令，但它似乎是消除数据依赖性的有益权衡。我假设寄存器越大增益越大（可能取决于芯片）。
更新 3
通过向量化 transition mask = ((~lt + gt) & lt) | ((~gt + lt) & gt)
我可以编译这个输出
vmovdqu ymm5,ymmword ptr transition_mask[rax]
vmovdqu ymm4,ymm5
vpandn ymm0,ymm5,ymm6
vpaddb ymm1,ymm0,ymm5
vpand ymm3,ymm1,ymm5
vpandn ymm2,ymm5,ymm6
vpaddb ymm0,ymm2,ymm5
vpand ymm1,ymm0,ymm5
vpor ymm3,ymm1,ymm3
vmovdqu ymmword ptr transition_mask[rax],ymm3
乍一看，与 post 处理（垂直扫描 + 附加到输出）的潜在条件相关陷阱相比，它似乎更有效，尽管处理 2 个流而不是 1 个流似乎是正确且合乎逻辑的.
这缺少每个周期生成初始状态的能力（从 0 到 1 或 -1 的转换）。
不确定是否有一种方法可以增强 transition_mask 生成“位旋转”，或者使用 auto initial _tzcnt_u32(mask0) > _tzcnt_u32(mask1)，就像 Soons 在此处使用的那样：
结论
@aqrit 分享的方法是使用改进的 bit-twiddling 解决方案为每个块加载找到转换，结果证明是运行时性能最高的方法。 热内循环只有 9 条 asm 指令长（每 2 个找到的项目与其他方法进行比较）使用 tzcnt 和 blsr 像这样
tzcnt rax,rcx
mov qword ptr [rbx+rdx*8],rax
blsr rcx,rcx
tzcnt rax,rcx
mov qword ptr [rbx+rdx*8+8],rax
blsr rcx,rcx
add rdx,2
cmp rdx,r8
jl main+2580h (...)
对于您的情况，完全矢量化不是最佳选择。这在技术上是可行的，但我认为生成 uint64_t 值数组的开销（我假设你正在为 64 位 CPU 编译）会吃掉所有的利润。
相反，您应该加载 32 字节的块，并立即将它们转换为位掩码。方法如下：
inline void loadBits( const int8_t* rsi, uint32_t& lt, uint32_t& gt )
{
const __m256i vec = _mm256_loadu_si256( ( const __m256i* )rsi );
lt = (uint32_t)_mm256_movemask_epi8( vec );
const __m256i cmp = _mm256_cmpgt_epi8( vec, _mm256_setzero_si256() );
gt = (uint32_t)_mm256_movemask_epi8( cmp );
}
您的其余代码应处理这些位图。要找到第一个 non-zero 元素（您只需在数据的开头这样做），扫描 (lt | gt) 整数中的最低有效设置位。要查找 -1 数字，扫描 lt 整数中的最低有效设置位，查找 +1 数字扫描 gt 整数中的最低有效设置位。找到并处理后，您可以使用按位与清除两个整数的低位部分，或者将它们都向右移动。
CPU 有 BSF instruction which scans for the lowest set bit in an integer, and returns two things at once: a flag telling if the integer was zero, and the index of that set bit. If you’re using VC++ there’s _BitScanForward 内在，否则使用内联 ASM，该指令仅在 VC++ 中公开； GCC 的 __builtin_ctz 不是一回事，它只是 returns 一个值而不是两个值。
但是，在 AMD CPU 上，TZCNT instruction from BMI 1 set is somewhat faster than the old-school BSF（在 Intel 上它们是相等的）。在 AMD 上，TZCNT 可能会稍微快一些，尽管有额外的指令与 0 进行比较。
在 64 位 SIMD 通道之间串行传输状态比在 64 位通用寄存器 (gpr) 之间串行传输状态更昂贵。
实际上，lookup-tables（或 SIMD left-packing）一次只能处理 8 个元素。如果数据平均每 64 个元素保留 6 个元素，那么 left-packing 会浪费大量处理 （特别是如果我们正在收集偏移量而不是进行收集操作）。如果位集密集，则考虑移动到 lookup-tables.
正如@Snoots 所建议的那样，使用 SIMD 创建 64 位位集并使用位扫描内在函数查找所需位集的索引。
分支预测错误：
使用 transition_mask = ((~lt + gt) & lt) | ((~gt + lt) & gt) 或@FalkHüffner transition_mask = (lt ^ (lt - gt)) & (gt ^ (gt – lt)).
其中一项算术运算的状态为 carry-in/carry-out。我会小心使用 _subborrow_u64，因为它是相当不常见的内在函数（并且在旧编译器上存在错误）。
剩下唯一的分支在位扫描操作上循环。必须提取所有设置位.. 但我们可以展开操作和超调以使分支更可预测。超调量需要调整到预期的数据集。
未测试。未检查 Asm。
#include <immintrin.h>
#include <stdint.h>
static inline
uint64_t get_mask (int8_t* src, unsigned char* state) {
__m256i src0 = _mm256_loadu_si256((__m256i*)(void*)src);
__m256i src1 = _mm256_loadu_si256((__m256i*)(void*)&src[32]);
uint64_t lt = (uint32_t)_mm256_movemask_epi8(src0) |
(((uint64_t)(uint32_t)_mm256_movemask_epi8(src1)) << 32);
src0 = _mm256_cmpgt_epi8(src0, _mm256_setzero_si256());
src1 = _mm256_cmpgt_epi8(src1, _mm256_setzero_si256());
uint64_t gt = (uint32_t)_mm256_movemask_epi8(src0) |
(((uint64_t)(uint32_t)_mm256_movemask_epi8(src1)) << 32);
// if borrow then greater-than span extends past the msb
uint64_t m;
unsigned char s = *state;
*state = _subborrow_u64(s, lt, gt, (unsigned long long*)&m); // sbb
return (m ^ lt) & ((gt - (lt + !s)) ^ gt);
}
static inline
size_t bitset_to_index (uint64_t* dst, uint64_t base, uint64_t mask) {
int64_t cnt = _mm_popcnt_u64(mask);
int64_t i = 0;
do { // unroll to taste...
dst[i + 0] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 1] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 2] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
dst[i + 3] = base + _tzcnt_u64(mask); mask = _blsr_u64(mask);
i += 4;
} while (i < cnt);
return (size_t)cnt;
}
static
uint64_t* get_transition_indices (uint64_t* dst, int8_t* src, size_t len) {
unsigned char state = 0; // in less-than span
uint64_t base = 0; // offset into src array
size_t end = len / 64;
for (size_t i = 0; i < end; i++) {
uint64_t mask = get_mask(src, &state);
src += 64;
dst += bitset_to_index(dst, base, mask);
base += 64;
}
if (len % 64) {
; // todo: tail loop
}
return dst;
}