在 SSE 寄存器中查找最常出现的元素
Posted
技术标签:
【中文标题】在 SSE 寄存器中查找最常出现的元素【英文标题】：Finding the most frequently occurring element in an SSE register 【发布时间】：2017-08-03 05:56:35 【问题描述】：有人对如何在 SSE4.x 中计算 8 位整数向量的模式（统计量）有任何想法吗？澄清一下，这将是 128 位寄存器中的 16x8 位值。
我希望将结果作为选择模式值元素的矢量掩码。即_mm_cmpeq_epi8(v, set1(mode(v)))的结果，以及标量值。
提供一些额外的上下文；虽然上述问题本身就是一个有趣的问题，但我已经通过了大多数我能想到的具有线性复杂度的算法。这个课程将抹去我计算这个数字所能获得的任何收益。
我希望在这里与大家一起寻找一些深奥的魔法。可能需要一个近似值来打破这个界限，例如“选择 a 频繁出现的元素”（NB 与 最 的差异），这将是优点。也可以使用概率答案。
SSE 和 x86 有一些非常有趣的语义。可能值得探索超优化通道。
【问题讨论】：
要找到模式，需要先计算直方图。但是由于内存冲突和依赖关系，很难有效地向量化直方图计算。可能的方法：***.com/questions/12985949 嗯...也许您可以通过将元素广播到新向量的每个成员来计算每个元素出现的频率，然后将pcmp 与原始向量和水平总和一起获得元素计数.现在使用pmaxub 找到最大值。不过好像也不算太快。
@fuz 幸运的是，水平总和不是必需的。这里可以使用更有效的垂直和。请参阅下面的答案。
@Martin Källman 有趣的问题。请注意，问题的数学标题较少可能有助于为您的问题吸引更多关注。例如：“如何有效地计算最常出现在 16x8 位 SSE 寄存器中的 8 位整数”。
@PeterCordes 在我的具体情况下，我将使用它为模式值元素创建掩码，用于其他目的
【参考方案1】：
这里可能适合使用相对简单的蛮力 SSEx 方法，请参见下面的代码。
想法是将输入向量 v 字节旋转 1 到 15 个位置并比较旋转后的向量
与原始的v 相等。缩短依赖链并增加
指令级并行性，两个计数器用于对这些相等的元素进行计数（垂直总和）：
sum1 和 sum2，因为可能有一些架构可以从中受益。
相等的元素计为 -1。变量sum = sum1 + sum2 包含带有值的总数
在 -1 和 -16 之间。 min_brc 包含广播到所有元素的sum 的水平最小值。
mask = _mm_cmpeq_epi8(sum,min_brc) 是被请求为模式值元素的掩码
OP的中间结果。在接下来的几行代码中，实际模式被提取出来。
此解决方案肯定比标量解决方案更快。 请注意，对于 AVX2，高 128 位通道可用于进一步加快计算速度。
仅计算模式值元素的掩码需要 20 个周期（吞吐量）。以实战模式播出 整个 SSE 寄存器大约需要 21.4 个周期。
注意下一个示例中的行为：
[1, 1, 3, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] 返回mask=[-1,-1,-1,-1,0,0,...,0]
并且众数值为 1，虽然 1 的出现频率与 3 一样。
下面的代码经过测试，但没有经过彻底测试
#include <stdio.h>
#include <x86intrin.h>
/* gcc -O3 -Wall -m64 -march=nehalem mode_uint8.c */
int print_vec_char(__m128i x);
__m128i mode_statistic(__m128i v)
__m128i sum2 = _mm_set1_epi8(-1); /* Each integer occurs at least one time */
__m128i v_rot1 = _mm_alignr_epi8(v,v,1);
__m128i v_rot2 = _mm_alignr_epi8(v,v,2);
__m128i sum1 = _mm_cmpeq_epi8(v,v_rot1);
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot2));
__m128i v_rot3 = _mm_alignr_epi8(v,v,3);
__m128i v_rot4 = _mm_alignr_epi8(v,v,4);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot3));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot4));
__m128i v_rot5 = _mm_alignr_epi8(v,v,5);
__m128i v_rot6 = _mm_alignr_epi8(v,v,6);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot5));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot6));
__m128i v_rot7 = _mm_alignr_epi8(v,v,7);
__m128i v_rot8 = _mm_alignr_epi8(v,v,8);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot7));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot8));
__m128i v_rot9 = _mm_alignr_epi8(v,v,9);
__m128i v_rot10 = _mm_alignr_epi8(v,v,10);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot9));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot10));
__m128i v_rot11 = _mm_alignr_epi8(v,v,11);
__m128i v_rot12 = _mm_alignr_epi8(v,v,12);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot11));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot12));
__m128i v_rot13 = _mm_alignr_epi8(v,v,13);
__m128i v_rot14 = _mm_alignr_epi8(v,v,14);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot13));
sum2 = _mm_add_epi8(sum2,_mm_cmpeq_epi8(v,v_rot14));
__m128i v_rot15 = _mm_alignr_epi8(v,v,15);
sum1 = _mm_add_epi8(sum1,_mm_cmpeq_epi8(v,v_rot15));
__m128i sum = _mm_add_epi8(sum1,sum2); /* Sum contains values such as -1, -2 ,...,-16 */
/* The next three instructions compute the horizontal minimum of sum */
__m128i sum_shft = _mm_srli_epi16(sum,8); /* Shift right 8 bits, while shifting in zeros */
__m128i min1 = _mm_min_epu8(sum,sum_shft); /* sum and sum_shuft are considered as unsigned integers. sum_shft is zero at the odd positions and so is min1 */
__m128i min2 = _mm_minpos_epu16(min1); /* Byte 0 within min2 contains the horizontal minimum of sum */
__m128i min_brc = _mm_shuffle_epi8(min2,_mm_setzero_si128()); /* Broadcast horizontal minimum */
__m128i mask = _mm_cmpeq_epi8(sum,min_brc); /* Mask = -1 at the byte positions where the value of v is equal to the mode of v */
/* comment next 4 lines out if there is no need to broadcast the mode value */
int bitmask = _mm_movemask_epi8(mask);
int indx = __builtin_ctz(bitmask); /* Index of mode */
__m128i v_indx = _mm_set1_epi8(indx); /* Broadcast indx */
__m128i answer = _mm_shuffle_epi8(v,v_indx); /* Broadcast mode to each element of answer */
/* Uncomment lines below to print intermediate results, to see how it works. */
// printf("sum = ");print_vec_char (sum );
// printf("sum_shft = ");print_vec_char (sum_shft );
// printf("min1 = ");print_vec_char (min1 );
// printf("min2 = ");print_vec_char (min2 );
// printf("min_brc = ");print_vec_char (min_brc );
// printf("mask = ");print_vec_char (mask );
// printf("v_indx = ");print_vec_char (v_indx );
// printf("answer = ");print_vec_char (answer );
return answer; /* or return mask, or return both .... :) */
int main()
/* To test throughput set throughput_test to 1, otherwise 0 */
/* Use e.g. perf stat -d ./a.out to test throughput */
#define throughput_test 0
/* Different test vectors */
int i;
char x1[16] = 5, 2, 2, 7, 21, 4, 7, 7, 3, 9, 2, 5, 4, 3, 5, 5;
char x2[16] = 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5;
char x3[16] = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;
char x4[16] = 1, 2, 3, 2, 1, 6, 7, 8, 2, 2, 2, 3, 3, 2, 15, 16;
char x5[16] = 1, 1, 3, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;
printf("\n15...0 = 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0\n\n");
__m128i x_vec = _mm_loadu_si128((__m128i*)x1);
printf("x_vec = ");print_vec_char(x_vec );
__m128i y = mode_statistic (x_vec);
printf("answer = ");print_vec_char(y );
#if throughput_test == 1
__m128i x_vec1 = _mm_loadu_si128((__m128i*)x1);
__m128i x_vec2 = _mm_loadu_si128((__m128i*)x2);
__m128i x_vec3 = _mm_loadu_si128((__m128i*)x3);
__m128i x_vec4 = _mm_loadu_si128((__m128i*)x4);
__m128i x_vec5 = _mm_loadu_si128((__m128i*)x5);
__m128i y1, y2, y3, y4, y5;
__asm__ __volatile__ ( "vzeroupper" : : : ); /* Remove this line on non-AVX processors */
for (i=0;i<100000000;i++)
y1 = mode_statistic (x_vec1);
y2 = mode_statistic (x_vec2);
y3 = mode_statistic (x_vec3);
y4 = mode_statistic (x_vec4);
y5 = mode_statistic (x_vec5);
x_vec1 = mode_statistic (y1 );
x_vec2 = mode_statistic (y2 );
x_vec3 = mode_statistic (y3 );
x_vec4 = mode_statistic (y4 );
x_vec5 = mode_statistic (y5 );
printf("mask mode = ");print_vec_char(y1 );
printf("mask mode = ");print_vec_char(y2 );
printf("mask mode = ");print_vec_char(y3 );
printf("mask mode = ");print_vec_char(y4 );
printf("mask mode = ");print_vec_char(y5 );
#endif
return 0;
int print_vec_char(__m128i x)
char v[16];
_mm_storeu_si128((__m128i *)v,x);
printf("%3hhi %3hhi %3hhi %3hhi | %3hhi %3hhi %3hhi %3hhi | %3hhi %3hhi %3hhi %3hhi | %3hhi %3hhi %3hhi %3hhi\n",
v[15],v[14],v[13],v[12],v[11],v[10],v[9],v[8],v[7],v[6],v[5],v[4],v[3],v[2],v[1],v[0]);
return 0;
输出：
15...0 = 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0
x_vec = 5 5 3 4 | 5 2 9 3 | 7 7 4 21 | 7 2 2 5
sum = -4 -4 -2 -2 | -4 -3 -1 -2 | -3 -3 -2 -1 | -3 -3 -3 -4
min_brc = -4 -4 -4 -4 | -4 -4 -4 -4 | -4 -4 -4 -4 | -4 -4 -4 -4
mask = -1 -1 0 0 | -1 0 0 0 | 0 0 0 0 | 0 0 0 -1
answer = 5 5 5 5 | 5 5 5 5 | 5 5 5 5 | 5 5 5 5
使用 Evgeny Kluev 的method 计算水平最小值。
【讨论】：
如果您从 set1(255) 开始，并为每场比赛添加-1 比较结果，我认为您可以在末尾使用 phminposuw。它对单词进行操作，因此您需要将pmovzxbw 低半部分和punpckhbw 高半部分用零。然后phminposuw 将为您提供每一半的最低计数器的索引。错误，那么你遇到了决定拿哪一半的问题。 pminuw 并使用pcmpeqw 找出哪个向量的计数值已经是最小值，然后使用该比较结果混合index_low 和index_high+8。这将是您原始向量的索引。
这并没有我想象的那么高效；它可能不会比你在 printf 之后的序列好多少，我认为它正在解决同样的问题。
您实际上是从 0 开始计数器，这对于您的结束代码来说很好，但 phminposuw 认为 0 小于 -3 = unsigned 253。您跳过将每个元素与其自身进行比较，这将使单例元素的计数为 -1。但当然，pcmpeqb same,same 是公认的用于生成全一的打破依赖关系的惯用语。 :P 不管怎样，phminposuw 清理增加了 2 条指令，以将计数器初始化为全 1。
@PeterCordes 确实，_mm_minpos_epu16 比我最初想象的更适合 8 位整数。另见here。我得想想。请注意，上面代码中的min4_brc = _mm_shuffle_epi8(min4,_mm_setzero_si128()) 是多余的，因为最小值已经在 16 个元素中。
@PeterCordes 确实phminposuw 比 shift/min 序列快。每个函数调用可以节省大约 2 个周期。请注意，吞吐量约为 21.4 个周期。对于 18 端口 5 操作，这看起来是合理的。 Skylake 上每个周期 3.71 条指令。【参考方案2】：
对寄存器中的数据进行排序。 插入排序可以在 16 (15) 个步骤中完成，通过将寄存器初始化为“Infinity”，它试图说明一个单调递减的数组并将新元素并行插入到所有可能的位置：
// e.g. FF FF FF FF FF FF FF FF FF FF FF FF FF FF FF 78
__m128i sorted = _mm_or_si128(my_array, const_FFFFF00);
for (int i = 1; i < 16; ++i)
// Trying to insert e.g. A0, we must shift all the FF's to left
// e.g. FF FF FF FF FF FF FF FF FF FF FF FF FF FF 78 00
__m128i shifted = _mm_bslli_si128(sorted, 1);
// Taking the MAX of shifted and 'A0 on all places'
// e.g. FF FF FF FF FF FF FF FF FF FF FF FF FF FF A0 A0
shifted = _mm_max_epu8(shifted, _mm_set1_epi8(my_array[i]));
// and minimum of the shifted + original --
// e.g. FF FF FF FF FF FF FF FF FF FF FF FF FF FF A0 78
sorted = _mm_min_epu8(sorted, shifted);
然后计算vec[n+1] == vec[n] 的掩码，将掩码移动到 GPR 并使用它来索引 32768 条目 LUT 以获得最佳索引位置。
在实际情况下，人们可能想要排序的不仅仅是一个向量；即一次对 16 个 16 项向量进行排序；
__m128i input[16]; // not 1, but 16 vectors
transpose16x16(input); // inplace vector transpose
sort(transpose); // 60-stage network exists for 16 inputs
// linear search -- result in 'mode'
__m128i mode = input[0];
__m128i previous = mode;
__m128i count = _mm_set_epi8(0);
__m128i max_count = _mm_setzero_si128(0);
for (int i = 1; i < 16; i++)
__m128i ¤t = input[i];
// histogram count is off by one
// if (current == previous) count++;
// else count = 0;
// if (count > max_count)
// mode = current, max_count = count
prev = _mm_cmpeq_epi8(prev, current);
count = _mm_and_si128(_mm_sub_epi8(count, prev), prev);
__m128i max_so_far = _mm_cmplt_epi8(max_count, count);
mode = _mm_blendv_epi8(mode, current, max_so_far);
max_count = _mm_max_epi8(max_count, count);
previous = current;
内部循环每个结果的总摊销成本为 7-8 条指令； 排序通常每个阶段有 2 条指令——即每个结果有 8 条指令，当 16 个结果需要 60 个阶段或 120 条指令时。 （这仍然将转置作为练习——但我认为它应该比排序快得多？）
因此，这应该在每个 8 位结果 24 条指令的范围内。
【讨论】：
通过 7 个插入阶段（同时到低 8 和高 8 值），然后是双音合并 8、4、2 和 1，可能会更好地平衡操作...... 这是一个非常简洁的插入排序 :) 排序是我想到的第一件事，但不幸的是，它的线性复杂性将抵消我在特定情况下的任何收益:( 我认为 16 个向量的字节元素转置会很昂贵（与排序相比并不简单）。只有 16 个寄存器（直到 AVX512），您将不得不溢出一些。除此之外，每个向量有 16 个元素需要做很多工作。shufps 是一个不错的 2 输入 shuffle，但仅适用于 4 个元素的块。 pblendvb 是 2 微指令，只能在 Intel HSW/BDW 的 port5 上运行。 （在其他 uarches 上没有那么糟糕；Skylake 可以在 1 uop 中为任何端口运行非 AVX 版本（在xmm0 中混合控制），或在 2 uop 中为任何端口运行 AVX-128 版本。）【参考方案3】：
用于与标量代码进行性能比较。主要部分未矢量化，但在表格清除和 tmp 初始化时矢量化。 （对于 fx8150，每次 f() 调用 168 个周期（22M 调用在 3.7 GHz 时在 1.0002 秒内完成））
#include <x86intrin.h>
unsigned char tmp[16]; // extracted values are here (single instruction, store_ps)
unsigned char table[256]; // counter table containing zeroes
char f(__m128i values)
_mm_store_si128((__m128i *)tmp,values);
int maxOccurence=0;
int currentValue=0;
for(int i=0;i<16;i++)
unsigned char ind=tmp[i];
unsigned char t=table[ind];
t++;
if(t>maxOccurence)
maxOccurence=t;
currentValue=ind;
table[ind]=t;
for(int i=0;i<256;i++)
table[i]=0;
return currentValue;
g++ 6.3 输出：
f: # @f
movaps %xmm0, tmp(%rip)
movaps %xmm0, -24(%rsp)
xorl %r8d, %r8d
movq $-15, %rdx
movb -24(%rsp), %sil
xorl %eax, %eax
jmp .LBB0_1
.LBB0_2: # %._crit_edge
cmpl %r8d, %esi
cmovgel %esi, %r8d
movb tmp+16(%rdx), %sil
incq %rdx
.LBB0_1: # =>This Inner Loop Header: Depth=1
movzbl %sil, %edi
movb table(%rdi), %cl
incb %cl
movzbl %cl, %esi
cmpl %r8d, %esi
cmovgl %edi, %eax
movb %sil, table(%rdi)
testq %rdx, %rdx
jne .LBB0_2
xorps %xmm0, %xmm0
movaps %xmm0, table+240(%rip)
movaps %xmm0, table+224(%rip)
movaps %xmm0, table+208(%rip)
movaps %xmm0, table+192(%rip)
movaps %xmm0, table+176(%rip)
movaps %xmm0, table+160(%rip)
movaps %xmm0, table+144(%rip)
movaps %xmm0, table+128(%rip)
movaps %xmm0, table+112(%rip)
movaps %xmm0, table+96(%rip)
movaps %xmm0, table+80(%rip)
movaps %xmm0, table+64(%rip)
movaps %xmm0, table+48(%rip)
movaps %xmm0, table+32(%rip)
movaps %xmm0, table+16(%rip)
movaps %xmm0, table(%rip)
movsbl %al, %eax
ret
【讨论】：
即使在 AMD CPU 上 int 和 SIMD ALU 不竞争相同的执行端口（就像它们在 Intel 上一样），前端仍然是这个版本的一个重要瓶颈。如果 @wim 在 Skylake 上的每次呼叫吞吐量为 21.4c 与您在 AMD Piledriver(?) 上获得的吞吐量接近，则您需要将大约 8 个 SIMD 向量交织到一个标量函数中。这似乎是合理的：palignr 的吞吐量在 Skylake 和 Bulldozer 系列上是相同的每时钟 1，而paddb/pcmpeqb 在 2 个不同的端口上运行。但是 SIMD 版本仍然会占用 3/4 的前端带宽，除了清理。
无论如何，也许可以将@wim 的 SIMD 版本交错到这个内部循环中，所以每 1 个整数执行 16 个 SIMD。或者可能在交错之前将内部循环展开 2，并且每 1 个整数执行 8 个 SIMD。这只会在具有多个线程的 Steamroller 上发挥作用，其中一对中的每个核心（共享一个 SIMD 单元）每个时钟可以解码 4 条指令。在 Bulldozer/Piledriver 上，每个“集群”的两个线程在一个解码器上获得交替循环，因此几乎没有额外的前端带宽来保持整数 ALU 的馈送而不会使 SIMD ALU 饿死。
所以你的意思是算法复杂度太高而收益太少（如果有的话，只在amd中）？
由于前端的限制，CPU 的执行单元比它们可以同时使用的要多。但是无论您运行什么工作负载，该工作负载都有很多执行单元，通常足以使用大部分前端问题/重命名容量。另请参阅en.wikipedia.org/wiki/Dark_silicon。
Bulldozer/Piledriver 的收益可能几乎与 Haswell/Skylake 一样少。或者没有 AVX，可能有一些额外的 MOVDQA 指令会占用@wim 版本中剩余的前端带宽，而不需要 ALU 微指令。不过，Steamroller 实际上很有趣，因为在一对内核上有 2 个线程，您的前端带宽比共享 SIMD 执行单元所能处理的要多得多。以上是关于在 SSE 寄存器中查找最常出现的元素的主要内容，如果未能解决你的问题，请参考以下文章