c++ - 如何加快LUT查找的直方图？

Question

首先，我有一个数组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都会用一个新数组初始化，我需要将 数组(每个包含1000个整数)中的所有1000个数组散列到(每个包含1000个整数)到1000个数组(每个包含32个整数)代表每个数组中有多少个整数。

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[1000][b[e]]++;

我可以在0.00237秒内映射1000 * 1000值。还有什么其他方法可以加快代码执行速度吗？ (就像SIMD吗？)这段代码是我程序的瓶颈。

Answer 1

这本质上是一个直方图问题。您正在使用16位查找表映射16位值和5位值，但之后只是对LUT结果进行直方图绘制。有关直方图的更多信息，请参见下文。

首先，您可以使用最小的数据类型来增加LUT(以及原始数据)的密度。在x86上，将零或符号扩展的8位或16位数据加载到寄存器中的成本几乎与常规32位int加载(假定都命中高速缓存)的成本相同，而8位或16位存储区也和32位存储区一样便宜。

由于您的数据大小超过了L1 d缓存大小(对于所有最新的Intel设计，均为32kiB)，并且您以分散模式访问它，因此可以通过减少缓存大小获得很多。 (有关x86性能的更多信息，请参见x86标签Wiki，尤其是Agner Fog的资料)。

由于a在每个平面上的条目数少于65536，因此bin计数永远不会溢出16位计数器，因此bins也可以是uint16_t。

您的copy()没有任何意义。为什么要复制到b[32768]而不是让内部循环使用指向当前LUT的指针？您以只读方式使用它。您仍然要复制的唯一原因是，如果您无法将产生不同LUT的代码更改为首先生成int或uin8_t，则将其从int8_t复制到uint8_t。

这个版本利用了这些想法和一些直方图技巧，并编译为看起来不错的asm(Godbolt compiler explorer: 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个即可，而不是对单个直方图进行 slice 。展开因子甚至不必为2的幂。

这样一来就无需在结尾处的 bins[i][k] = sum(tmpbins[0..3][k])上进行 k循环。

在使用前将 bins[i..i+unroll_factor][0..31]归零，而不是在循环外将整个事物归零。这样，当您启动时，所有垃圾箱在L1缓存中都会很热，并且这项工作可能与内部循环中比较繁重的工作重叠。

硬件预取器可以跟踪多个顺序流，因此不必担心从 a加载时会有更多的缓存丢失。 (为此也要使用 vector 加载，并在加载后将它们切成薄片)。

其他有关直方图的有用答案的问题:

Methods to vectorise histogram in SIMD?建议使用多个bin数组并在最后的技巧上求和。

Optimizing SIMD histogram calculation x86 asm加载a值的 vector ，并使用pextrb提取到整数寄存器。 (在您的代码中，您将使用 pextrw / _mm_extract_epi16 )。 在所有加载/存储操作发生的情况下，进行 vector 加载并使用ALU ops进行解压缩很有意义。如果L1命中率很高，则内存uop吞吐量可能是瓶颈，而不是内存/缓存延迟。

How to optimize histogram statistics with neon intrinsics?有一些相同的想法:bins数组的多个副本。它还具有ARM特定的建议，用于在SIMD vector 中执行地址计算(ARM可以在一条指令中从 vector 中获得两个标量)，并以相反的方式布置多仓数组。

AVX2收集LUT的说明

如果要在Intel Skylake上运行此程序，甚至可以使用AVX2进行LUT查找以收集指令。 (在Broadwell上，这可能是一个收支平衡，而在Haswell上，它将失去；它们支持 vpgatherdd ( _mm_i32gather_epi32 )，但实现效率不高。希望Skylake避免在元素之间存在重叠时多次击中同一条缓存行)。

是的，即使最小的收集粒度是32位元素，您仍然可以从 uint16_t数组(比例因子= 2)进行收集。这意味着您会在每个32位 vector 元素的上半部分得到垃圾，而不是0，但这没关系。高速缓存行拆分并不理想，因为我们可能在高速缓存吞吐量方面遇到瓶颈。

收集到的元素的上半部分中的垃圾无关紧要，因为无论如何您都只使用 pextrw提取了有用的16位。 (并使用标量代码执行过程的直方图部分)。

只要每个元素都来自直方图箱的单独 slice /平面，就可以潜在地使用另一个聚集从直方图箱中加载。否则，如果两个元素来自同一bin，则只有当您手动将递增的 vector 散布回直方图(带有标量存储)时，它才会递增1。这种针对散布店的冲突检测是 AVX512CD存在的原因。 AVX512确实具有分散指令以及聚集指令(已在AVX2中添加)。

AVX512

请参见 page 50 of Kirill Yukhin's slides from 2014，以获取一个示例循环，该循环将重试直到没有冲突为止。但它没有显示如何用 get_conflict_free_subset()( __m512i _mm512_conflict_epi32 (__m512i a))(它与之冲突的所有先前元素的每个元素返回一个位图)来实现 vpconflictd。正如@Mysticial指出的那样，简单的实现要比冲突检测指令仅产生掩码寄存器结果而不是另一个 vector 要简单。

我进行了搜索，但没有找到英特尔发布的有关使用AVX512CD的教程/指南，但是大概他们认为在某些情况下在 _mm512_lzcnt_epi32的结果上使用 vplzcntd( vpconflictd)是有用的，因为它也是AVX512CD的一部分。

也许您“应该”做一些更聪明的事情，而不只是跳过所有有冲突的元素？也许可以检测到标量回退会更好的情况，例如所有16个dword元素都具有相同的索引？ vpbroadcastmw2d将掩码寄存器广播到结果的所有32位元素，因此您可以将mask-register值与 vpconflictd中每个元素中的位图对齐。 (并且AVX512F中的元素之间已经存在比较，按位和其他操作)。

Kirill的幻灯片列出了 VPTESTNM{D,Q}(来自AVX512F)以及冲突检测指令。它从 DEST[j] = (SRC1[i+31:i] BITWISE AND SRC2[i+31:i] == 0)? 1 : 0生成一个掩码。即将AND元素放在一起，如果不相交，则将该元素的 mask 结果设置为1。

可能也相关: http://colfaxresearch.com/knl-avx512/表示 "For a practical illustration, we construct and optimize a micro-kernel for particle binning particles"，带有一些AVX2的代码(我认为)。但这是我没有做的免费注册的背后。根据该图，我认为他们在将一些矢量化的东西生成想要分散的数据之后，将实际的分散部分作为标量进行处理。第一个链接表示第二个链接“用于先前的指令集”。