iou*_*vxz 4 c++ optimization simd histogram
首先，我有一个数组int a[1000][1000]。所有这些整数都在 0 到 32767 之间，它们是已知的常量：它们在程序运行期间永远不会改变。
其次，我有一个数组 b[32768]，它包含 0 到 32 之间的整数。我使用这个数组将 a 中的所有数组映射到 32 个 bin：
int bins[32]{};
for (auto e : a[i])//mapping a[i] to 32 bins.
bins[b[e]]++;
每次，数组 b 将用一个新数组初始化，我需要将数组 a 中的所有 1000 个数组（每个包含 1000 个整数）散列到 1000 个数组，每个数组包含 32 个整数，表示有多少整数落入其每个 bin 。
int new_array[32768] = {some new mapping};
copy(begin(new_array), end(new_array), begin(b));//reload array b;
int bins[1000][32]{};//output array to store results .
for (int i = 0; i < 1000;i++)
for (auto e : a[i])//hashing a[i] to 32 bins.
bins[i][b[e]]++;
我可以在 0.00237 秒内映射 1000*1000 个值。有没有其他方法可以加速我的代码？（比如SIMD？）这段代码是我程序的瓶颈。
这本质上是一个直方图问题。您使用 32k 条目的查找表将 16 位值映射到 5 位值，但之后它只是对 LUT 结果进行直方图绘制。像++counts[ b[a[j]] ];，这里counts是bins[i]。有关直方图的更多信息，请参见下文。
首先，您可以使用尽可能小的数据类型来增加 LUT（和原始数据）的密度。在 x86 上，将 8 位或 16 位数据加载到寄存器中的零或符号扩展加载与常规 32 位int加载（假设两者都命中缓存）的成本几乎完全相同，而 8 位或 16 -bit store 也和 32-bit store 一样便宜。
由于您的数据大小超过 L1 d-cache 大小（对于所有最新的 Intel 设计为 32kiB），并且您以分散的模式访问它，因此您可以从缩小缓存占用空间中获益良多。（有关 x86 性能的更多信息，请参阅x86标签 wiki，尤其是 Agner Fog 的内容）。
由于a每个平面中的条目少于 65536 个，因此您的 bin 计数永远不会溢出 16 位计数器，因此bins也可以uint16_t。
你的copy()毫无意义。为什么要复制b[32768]而不是让内部循环使用指向当前 LUT 的指针？您以只读方式使用它。你仍然要复制的唯一原因是从复制int到uin8_t如果你不能改变会产生不同的LUT，产生的代码int8_t或uint8_t摆在首位。
这个版本利用了这些想法和一些直方图技巧，并编译为看起来不错的 asm（Godbolt 编译器资源管理器：gcc6.2 -O3 -march=haswell(AVX2)）：
// untested
//#include <algorithm>
#include <stdint.h>
const int PLANES = 1000;
void use_bins(uint16_t bins[PLANES][32]); // pass the result to an extern function so it doesn't optimize away
// 65536 or higher triggers the static_assert
alignas(64) static uint16_t a[PLANES][1000]; // static/global, I guess?
void lut_and_histogram(uint8_t __restrict__ lut[32768])
{
alignas(16) uint16_t bins[PLANES][32]; // don't zero the whole thing up front: that would evict more data from cache than necessary
// Better would be zeroing the relevant plane of each bin right before using.
// you pay the rep stosq startup overhead more times, though.
for (int i = 0; i < PLANES;i++) {
alignas(16) uint16_t tmpbins[4][32] = {0};
constexpr int a_elems = sizeof(a[0])/sizeof(uint16_t);
static_assert(a_elems > 1, "someone changed a[] into a* and forgot to update this code");
static_assert(a_elems <= UINT16_MAX, "bins could overflow");
const uint16_t *ai = a[i];
for (int j = 0 ; j<a_elems ; j+=4) { //hashing a[i] to 32 bins.
// Unrolling to separate bin arrays reduces serial dependencies
// to avoid bottlenecks when the same bin is used repeatedly.
// This has to be balanced against using too much L1 cache for the bins.
// TODO: load a vector of data from ai[j] and unpack it with pextrw.
// even just loading a uint64_t and unpacking it to 4 uint16_t would help.
tmpbins[0][ lut[ai[j+0]] ]++;
tmpbins[1][ lut[ai[j+1]] ]++;
tmpbins[2][ lut[ai[j+2]] ]++;
tmpbins[3][ lut[ai[j+3]] ]++;
static_assert(a_elems % 4 == 0, "unroll factor doesn't divide a element count");
}
// TODO: do multiple a[i] in parallel instead of slicing up a single run.
for (int k = 0 ; k<32 ; k++) {
// gcc does auto-vectorize this with a short fully-unrolled VMOVDQA / VPADDW x3
bins[i][k] = tmpbins[0][k] + tmpbins[1][k] +
tmpbins[2][k] + tmpbins[3][k];
}
}
// do something with bins. An extern function stops it from optimizing away.
use_bins(bins);
}
内循环 asm 如下所示：
.L2:
movzx ecx, WORD PTR [rdx]
add rdx, 8 # pointer increment over ai[]
movzx ecx, BYTE PTR [rsi+rcx]
add WORD PTR [rbp-64272+rcx*2], 1 # memory-destination increment of a histogram element
movzx ecx, WORD PTR [rdx-6]
movzx ecx, BYTE PTR [rsi+rcx]
add WORD PTR [rbp-64208+rcx*2], 1
... repeated twice more
使用来自 rbp 的 32 位偏移量（而不是来自 rsp 的 8 位偏移量，或使用另一个寄存器：/），代码密度并不好。尽管如此，平均指令长度并不会太长，以至于它可能会成为任何现代 CPU 上指令解码的瓶颈。
由于无论如何您都需要制作多个直方图，因此只需并行制作 4 到 8 个直方图，而不是将 bin 切片以获得单个直方图。展开因子甚至不必是 2 的幂。
这消除了最后bins[i][k] = sum(tmpbins[0..3][k])循环的需要k。
零bins[i..i+unroll_factor][0..31]权的使用之前，而不是零，循环外的整个事情。这样一来，当您开始时，L1 缓存中的所有 bin 都会很热，并且这项工作可能会与内部循环的负载更重的工作重叠。
硬件预取器可以跟踪多个顺序流，因此不必担心从a. （为此也使用矢量加载，并在加载后将它们切碎）。
其他有关直方图有用答案的问题：
a值向量并使用 pextrb 提取到整数寄存器。（在您的代码中，您将使用pextrw/_mm_extract_epi16）。 随着所有加载/存储 uops 的发生，进行向量加载并使用 ALU 操作来解包是有意义的。具有良好的 L1 命中率，内存 uop 吞吐量可能是瓶颈，而不是内存/缓存延迟。如果您打算在 Intel Skylake 上运行它，您甚至可以使用 AVX2 收集指令进行 LUT 查找。（在 Broadwell 上，它可能是一个盈亏平衡点，而在 Haswell 上它会失败；它们支持vpgatherdd( _mm_i32gather_epi32)，但没有那么高效的实现。希望 Skylake 在元素之间存在重叠时避免多次访问相同的缓存行） .
是的，uint16_t即使最小的收集粒度是 32 位元素，您仍然可以从（比例因子 = 2）的数组中收集。这意味着您在每个 32 位向量元素的高半部分得到垃圾而不是 0，但这无关紧要。缓存行拆分并不理想，因为我们可能会遇到缓存吞吐量的瓶颈。
收集元素的前半部分中的垃圾无关紧要，因为无论如何您都只提取有用的 16 位pextrw。（并使用标量代码执行过程的直方图部分）。
您可以潜在地使用另一个收集从直方图箱加载，只要每个元素来自直方图箱的单独切片/平面。否则，如果两个元素来自同一个 bin，当您手动将增加的向量散布回直方图（使用标量存储）时，它只会增加 1。这种针对分散存储的冲突检测是AVX512CD存在的原因。AVX512 确实有分散指令，以及收集（已经在 AVX2 中添加）。
AVX512
请参阅Kirill Yukhin 2014 年幻灯片的第 50 页，了解重试直到没有冲突的示例循环；但它没有显示如何get_conflict_free_subset()在__m512i _mm512_conflict_epi32 (__m512i a)( vpconflictd)方面实现（它在与它冲突的所有先前元素的每个元素中返回一个位图）。正如@Mysticial 指出的那样，如果冲突检测指令只是产生一个掩码寄存器结果，而不是另一个向量，那么简单的实现就没有那么简单。
我搜索了但没有找到英特尔发布的关于使用 AVX512CD 的教程/指南，但大概他们认为在结果上使用_mm512_lzcnt_epi32( vplzcntd) 在vpconflictd某些情况下很有用，因为它也是 AVX512CD 的一部分。
也许你“应该”做一些比跳过所有有任何冲突的元素更聪明的事情？也许检测标量回退会更好的情况，例如所有 16 个双字元素都具有相同的索引？ vpbroadcastmw2d将掩码寄存器广播到结果的所有 32 位元素，这样您就可以将掩码寄存器值与vpconflictd. （并且已经有来自 AVX512F 的元素之间的比较、按位和其他操作）。
Kirill 的幻灯片列表VPTESTNM{D,Q}（来自 AVX512F）以及冲突检测说明。它从 生成一个掩码DEST[j] = (SRC1[i+31:i] BITWISE AND SRC2[i+31:i] == 0)? 1 : 0。即AND元素在一起，如果它们不相交，则将该元素的掩码结果设置为1。
可能也相关：http : //colfaxresearch.com/knl-avx512/说“为了实际说明，我们为粒子合并粒子构建和优化了一个微内核”，以及 AVX2 的一些代码（我认为）。但这背后是我没有做过的免费注册。根据图表，我认为他们正在将实际的散布部分作为标量进行处理，经过一些矢量化的东西来生成他们想要散布的数据。第一个链接说第二个链接是“用于以前的指令集”。
当存储桶的数量与数组的大小相比较小时，复制计数数组并展开以最小化重复元素的存储转发延迟瓶颈变得可行。但是对于聚集/分散策略，如果我们为每个向量元素使用不同的数组，它也避免了冲突的可能性，解决了正确性问题。
当收集/分散指令仅采用一个数组基数时，我们如何做到这一点？使所有计数数组连续，并用一条额外的 SIMD 添加指令偏移每个索引向量，完全取代冲突检测和分支。
如果桶的数量不是 16 的倍数，您可能想要四舍五入数组几何，以便每个计数子集从对齐的偏移量开始。或者不，如果缓存局部性比避免最终减少的未对齐负载更重要。
const size_t nb = 32; // number of buckets
const int VEC_WIDTH = 16; // sizeof(__m512i) / sizeof(uint32_t)
alignas(__m512i) uint32_t counts[nb * VEC_WIDTH] = {0};
// then in your histo loop
__m512i idx = ...; // in this case from LUT lookups
idx = _mm512_add_epi32(idx, _mm512_setr_epi32(
0*nb, 1*nb, 2*nb, 3*nb, 4*nb, 5*nb, 6*nb, 7*nb,
8*nb, 9*nb, 10*nb, 11*nb, 12*nb, 13*nb, 14*nb, 15*nb));
// note these are C array indexes, not byte offsets
__m512i vc = _mm512_i32gather_epi32(idx, counts, sizeof(counts[0]));
vc = _mm512_add_epi32(vc, _mm512_set1_epi32(1));
_mm512_i32scatter_epi32(counts, idx, vc, sizeof(counts[0]));
https://godbolt.org/z/8Kesx7sEK显示上面的代码实际上可以编译。（在循环内，向量常数设置可能会被提升，但不会在每次收集或分散之前将掩码寄存器设置为全一，或准备归零合并目标。）
然后在主直方图循环之后，减少到一个计数数组：
// Optionally with size_t nb as an arg
// also optionally use restrict if you never reduce in-place, into the bottom of the input.
void reduce_counts(int *output, const int *counts)
{
for (int i = 0 ; i < nb - (VEC_WIDTH-1) ; i+=VEC_WIDTH) {
__m512i v = _mm512_load_si512(&counts[i]); // aligned load, full cache line
// optional: unroll this and accumulate two vectors in parallel for better spatial locality and more ILP
for (int offset = nb; offset < nb*VEC_WIDTH ; offset+=nb) {
__m512i tmp = _mm512_loadu_si512(&counts[i + offset]);
v = _mm512_add_epi32(v, tmp);
}
_mm512_storeu_si512(&output[i], v);
}
// if nb isn't a multiple of the vector width, do some cleanup here
// possibly using a masked store to write into a final odd-sized destination
}
显然这对于太多的桶是不好的；你最终不得不将更多的内存归零，并在最后循环很多。使用 256 位而不是 512 位收集会有所帮助，您只需要一半的数组，但是收集/分散指令的效率会随着向量的增加而提高。例如vpgatherdd，在 Cascade Lake 上，256 位每 5 个周期一个，512 位每 9.25 个周期一个。（两者都是前端的 4 uop）
或者在 Ice Lake 上，vpscatterdd每 7 个周期一个 ymm，每 11 个周期一个 zmm。（与 2x ymm 的 14 相比）。 https://uops.info/
在你的bins[1000][32]情况下，你实际上可以使用后面的元素bins[i+0..15]作为额外的计数数组，如果你先清零，至少在前 1000-15 次外循环迭代中。这将避免触及额外的内存：下一个外部循环的归零将从前一个counts[32]有效开始。
（这对于 C 2D 和 1D 数组来说会有点快和松散，但是在[32]C 数组类型末尾之后的所有实际访问都将通过memset（即unsigned char*）或通过_mm*也允许别名任何东西的内在函数）
有关的：
count[0] += (arr[i] == 0)等等，您可以使用 SIMD 打包比较对其进行矢量化 -大型数组或列表的 4 桶直方图的微优化 当桶数小于或等于时这很有趣SIMD 向量中的元素数。