缓存友好的离线随机读取

Question

考虑一下 C++ 中的这个函数：

void foo(uint32_t *a1, uint32_t *a2, uint32_t *b1, uint32_t *b2, uint32_t *o) {
    while (b1 != b2) {
        // assert(0 <= *b1 && *b1 < a2 - a1)
        *o++ = a1[*b1++];
    }
}

它的目的应该足够明确。不幸的是，b1 包含 随机 数据并丢弃缓存，使foo 成为我程序的瓶颈。无论如何我可以优化它吗？

这是一个 SSCCE，应该类似于我的实际代码：

#include <iostream>
#include <chrono>
#include <algorithm>
#include <numeric>

namespace {
    void foo(uint32_t *a1, uint32_t *a2, uint32_t *b1, uint32_t *b2, uint32_t *o) {
        while (b1 != b2) {
            // assert(0 <= *b1 && *b1 < a2 - a1)
            *o++ = a1[*b1++];
        }
    }

    constexpr unsigned max_n = 1 << 24, max_q = 1 << 24;
    uint32_t data[max_n], index[max_q], result[max_q];
}

int main() {
    uint32_t seed = 0;
    auto rng = [&seed]() { return seed = seed * 9301 + 49297; };
    std::generate_n(data, max_n, rng);
    std::generate_n(index, max_q, [rng]() { return rng() % max_n; });

    auto t1 = std::chrono::high_resolution_clock::now();
    foo(data, data + max_n, index, index + max_q, result);
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << std::chrono::duration<double>(t2 - t1).count() << std::endl;

    uint32_t hash = 0;
    for (unsigned i = 0; i < max_q; i++)
        hash += result[i] ^ (i << 8) ^ i;
    std::cout << hash << std::endl;
}

这不是Cache-friendly copying of an array with readjustment by known index, gather, scatter，它询问随机写入并假设b 是一个排列。

Answer 1

首先我们来看看上面代码的实际表现：

$ sudo perf stat ./offline-read
0.123023
1451229184

 Performance counter stats for './offline-read':

        184.661547      task-clock (msec)         #    0.997 CPUs utilized          
                 3      context-switches          #    0.016 K/sec                  
                 0      cpu-migrations            #    0.000 K/sec                  
               717      page-faults               #    0.004 M/sec                  
       623,638,834      cycles                    #    3.377 GHz                    
       419,309,952      instructions              #    0.67  insn per cycle         
        70,803,672      branches                  #  383.424 M/sec                  
            16,895      branch-misses             #    0.02% of all branches        

       0.185129552 seconds time elapsed

我们得到了 0.67 的低 IPC，这可能几乎完全是由于 DRAM⁵ 的加载缺失造成的。让我们确认一下：

sudo ../pmu-tools/ocperf.py stat -e cycles,LLC-load-misses,cycle_activity.stalls_l3_miss ./offline-read
perf stat -e cycles,LLC-load-misses,cpu/event=0xa3,umask=0x6,cmask=6,name=cycle_activity_stalls_l3_miss/ ./offline-read
0.123979
1451229184

 Performance counter stats for './offline-read':

       622,661,371      cycles                                                      
        16,114,063      LLC-load-misses                                             
       368,395,404      cycle_activity_stalls_l3_miss                                   

       0.184045411 seconds time elapsed

因此，620k 中的约 370k 周期直接因未完成的未命中而停滞。事实上，foo() 中以这种方式停滞的周期部分要高得多，接近 90%，因为 perf 也在测量 init 和 accumulate 代码，这需要大约三分之一的运行时间（但没有重要的 L3 未命中）。

这并不意外，因为我们知道随机读取模式a1[*b1++] 的局部性基本上为零。其实LLC-load-misses的数量是1600万¹，几乎正好对应a1的1600万随机读取。²

如果我们假设 foo() 的 100% 都在等待内存访问，我们可以了解每次未命中的总成本：0.123 sec / 16,114,063 misses == 7.63 ns/miss。在我的机器上，内存延迟在最佳情况下约为 60 ns，因此每次未命中少于 8 ns 意味着我们已经在提取大量内存级并行 (MLP)：大约 8 个未命中必须重叠并且在-平均飞行来实现这一点（甚至完全忽略来自b1 的流式负载和o 的流式写入的额外流量）。

所以我不认为有很多调整可以应用于简单循环来做得更好。不过，有两种可能性：

• 写入o 的非临时存储（如果您的平台支持）。这将消除RFO 对普通商店所暗示的读取。这应该是一个直接的胜利，因为o 再也不会被读取（在计时部分内！）。

• 软件预取。仔细调整 a1 或 b1 的预取可能会有所帮助有点。然而，影响将是相当有限的，因为我们已经接近上述 MLP 的极限。此外，我们预计b1 的线性读取几乎可以被硬件预取器完美地预取。 a1 的随机读取似乎可以进行预取，但实际上循环中的 ILP 通过乱序处理（至少在像最近的 x86 这样的大型 OoO 处理器上）会导致足够的 MLP。

  在 cmets 用户中，harold 已经提到他尝试预取，但效果很小。

因此，由于简单的调整不太可能产生太大的效果，因此您只需要转换循环即可。一种“明显”的转换是对索引b1（连同索引元素的原始位置）进行排序，然后按排序顺序从a1 进行读取。这将a1 的读取从完全随机转换为几乎³ 线性，但现在写入都是随机的，这也好不到哪里去。

排序然后取消排序

关键问题是a1 在b1 控制下的读取是随机的，而a1 很大，基本上每次读取都会导致 DRAM 未命中。我们可以通过对b1 进行排序来解决这个问题，然后读取a1 以获得置换结果。现在您需要“取消排列”结果 a1 以获得最终顺序的结果，这只是另一种排序，这次是在“输出索引”上。

这是一个工作示例，给定输入数组 a、索引数组 b 和输出数组 o，以及 i，这是每个元素的（隐式）位置：

  i =   0   1   2   3
  a = [00, 10, 20, 30]
  b = [ 3,  1,  0,  1]
  o = [30, 10, 00, 10] (desired result)

首先，对数组b进行排序，将原始数组位置i作为辅助数据（或者您可以将其视为排序元组(b[0], 0), (b[1], 1), ...），这将为您提供排序后的b数组b'和排序后的索引列表i'如图：

  i' = [ 2, 1, 3, 0]
  b' = [ 0, 1, 1, 3]

现在您可以在b' 的控制下从a 读取置换结果数组o'。此读取按顺序严格增加，并且应该能够以接近memcpy 的速度运行。事实上，您可以利用广泛的连续 SIMD 读取和一些 shuffle 来执行多次读取，并且一次将 4 字节元素移动到正确的位置（复制一些元素并跳过其他元素）：

  a  = [00, 10, 20, 30]
  b' = [ 0,  1,  1,  3]
  o' = [00, 10, 10, 30]

最后，您对o' 进行反置换以得到o，从概念上讲，只需在置换索引i' 上对o' 进行排序即可：

  i' = [ 2,  1,  3,  0]
  o' = [00, 10, 10, 30]
  i  = [ 0,  1,  2,  3]
  o  = [30, 10, 00, 10]

完成！

现在这是该技术中最简单的想法，并且对缓存不是特别友好（每次传递概念上迭代一个或多个 2^26 字节数组），但它至少完全使用它读取的每个缓存行（不像原始循环仅从缓存行中读取单个元素，这就是为什么即使数据仅占用 100 万个缓存行也会有 1600 万次未命中的原因！）。所有的读取或多或少都是线性的，所以硬件预取会有很大帮助。

您可能获得多大的加速取决于您将如何实现排序：它们需要快速且缓存敏感。几乎可以肯定，某种类型的缓存感知基数排序效果最好。

以下是一些关于如何进一步改进的说明：

优化排序数量

您实际上不需要对b 进行完全排序。您只想对它进行“足够”排序，以便在b' 的控制下对a 的后续读取或多或少是线性的。例如，16 个元素适合一个缓存行，因此您根本不需要根据最后 4 位进行排序：无论如何都会读取相同的缓存行线性序列。您还可以对更少的位进行排序：例如，如果您忽略了 5 个最低有效位，您将以“几乎线性”的方式读取缓存行，有时会从完美线性模式中交换两条缓存行，例如：0, 1, 3, 2, 5, 4, 6, 7 .在这里，您仍然可以获得 L1 缓存的全部好处（对缓存行的后续读取总是会命中），我怀疑这种模式仍然可以很好地预取，如果不是，您总是可以通过软件预取来帮助它。

您可以在系统上测试最佳忽略位数是多少。忽略位有两个好处：

• 在基数搜索中要做的工作更少，无论是通过更少的通道，还是在一个或多个通道中需要更少的桶（这有助于缓存）。
• 在最后一步中“撤消”排列可能要做的工作更少：如果通过检查原始索引数组 b 进行撤消，则忽略位意味着您在撤消搜索时获得相同的节省。

缓存块的工作

上面的描述在几个连续的、不相交的通道中列出了所有内容，每个通道都在整个数据集上工作。在实践中，您可能希望将它们交错以获得更好的缓存行为。例如，假设您使用 MSD radix-256 排序，您可能会执行第一遍，将数据排序到 256 个桶中，每个桶约有 256K 元素。

然后，您可以只完成前几个（或前几个）存储桶的排序，而不是执行完整的第二遍，然后根据生成的b' 块继续读取a。您可以保证此块是连续的（即最终排序序列的后缀），因此您不会放弃读取中的任何位置，并且您的读取通常会被缓存。您也可以执行第一次去置换o'，因为o' 的块在缓存中也很热（也许您可以将后两个阶段组合成一个循环）。

智能反排列

优化的一个方面是如何精确地实现o' 的去排列。在上面的描述中，我们假设一些索引数组i 最初的值是[0, 1, 2, ..., max_q]，它与b 一起排序。这就是概念上的工作原理，但您可能不需要立即实际具体化i 并将其作为辅助数据进行排序。例如，在基数排序的第一遍中，i 的值是隐式已知的（因为您正在遍历数据），因此可以免费计算它⁴ 并在第一遍没有每个都按排序顺序出现。

可能还有比维护完整索引更有效的方法来执行“取消排序”操作。例如，原始未排序的b 数组在概念上具有执行 unsort 所需的所有信息，但我很清楚如何使用它来有效地 unsort。

会不会更快？

那么这实际上会比天真的方法更快吗？它在很大程度上取决于实现细节，尤其是实现排序的效率。在我的硬件上，简单的方法是每秒处理大约 1.4 亿个元素。缓存感知基数排序的在线描述似乎从每秒 200 到 6 亿个元素不等，并且由于您需要其中两个，如果您相信这些数字，那么大幅加速的机会似乎有限。另一方面，这些数字来自较旧的硬件，用于更一般的搜索（例如，对于所有 32 位密钥，而我们可能只使用 16 位）。

只有仔细实施才能确定是否可行，而可行性也取决于硬件。例如，在不能支持那么多 MLP 的硬件上，sorting-unsorting 方法变得相对更有利。

最佳方法还取决于max_n 和max_q 的相对值。例如，如果max_n >> max_q，那么即使使用最佳排序，读取也将是“稀疏”的，因此天真的方法会更好。另一方面，如果max_n << max_q，则通常会多次读取相同的索引，因此排序方法将具有良好的读取局部性，排序步骤本身将具有更好的局部性，并且可以进一步优化显式处理重复读取.

多核

您是否有兴趣并行化这个问题并不清楚。 foo() 的简单解决方案已经承认“直接”并行化，您只需将 a 和 b 数组划分为相等大小的块，针对每个线程，这似乎提供了完美的加速。不幸的是，您可能会发现您比线性扩展更糟糕，因为您将遇到内存控制器中的资源争用和关联的非核心/非核心资源，这些资源在套接字上的所有核心之间共享。因此，当您添加更多内核⁶时，对于内存的纯并行随机读取负载，您将获得多少吞吐量并不清楚。

对于基数排序版本，大多数瓶颈（存储吞吐量、总指令吞吐量）都在内核中，因此我希望它可以通过增加内核来合理扩展。正如 Peter 在评论中提到的，如果您使用超线程，则排序可能具有核心本地 L1 和 L2 缓存中良好局部性的额外好处，有效地让每个兄弟线程使用整个缓存，而不是将有效容量减少一半.当然，这涉及到仔细管理您的线程关联性，以便同级线程实际使用附近的数据，而不仅仅是让调度程序做任何事情。

---

¹ 你可能会问为什么 LLC-load-misses 不是说 32 或 4800 万，因为我们还必须读取 b1 的所有 1600 万个元素，然后是 accumulate() 调用读取所有result。答案是LLC-load-misses 仅计算在 L3 中实际上未命中的需求未命中。其他提到的读取模式是完全线性的，因此预取器将始终在需要之前将行带入 L3。根据 perf 的定义，这些不算作“LLC 未命中”。

²你可能想知道我是如何知道负载缺失都来自foo中a1的读取：我只是使用了perf record和perf mem 确认未命中来自预期的汇编指令。

³几乎线性，因为b1不是所有索引的排列，所以原则上可以跳过和重复索引。然而，在高速缓存行级别，很有可能每个高速缓存行都将按顺序读取，因为每个元素都有约 63% 的机会被包含在内，并且一个高速缓存行有 16 个 4 字节元素，所以只有任何给定的缓存都有零个元素的可能性大约为千万分之一。因此，在缓存行级别工作的预取可以正常工作。

⁴ 这里我的意思是价值的计算是免费的或者几乎是免费的，但当然写入仍然是收费的。然而，这仍然比“预先实现”方法好得多，它首先创建 i 数组 [0, 1, 2, ...] 需要 max_q 写入，然后再次需要另一个 max_q 写入以在第一次基数排序中对其进行排序经过。隐式物化只会导致第二次写入。

⁵ 事实上，实际定时部分foo() 的IPC 比低很多：根据我的计算，大约为0.15。整个过程上报的IPC是定时段的IPC和前后初始化累加代码的平均，IPC要高很多。

⁶ 值得注意的是，这与依赖负载延迟绑定工作流的扩展方式不同：负载正在执行随机读取但只能进行一个负载，因为每个负载都取决于结果的最后一个可以很好地扩展到多个内核，因为负载的串行性质不使用很多下游资源（但是从概念上讲，即使在单个内核上也可以通过更改内核循环来处理多个相关的负载流，从而加快此类负载的速度并行）。

Answer 2

您可以将索引分区到索引的较高位相同的存储桶中。请注意，如果索引不是随机的，则存储桶将溢出。

#include <iostream>
#include <chrono>
#include <cassert>
#include <algorithm>
#include <numeric>
#include <vector>

namespace {
constexpr unsigned max_n = 1 << 24, max_q = 1 << 24;

void foo(uint32_t *a1, uint32_t *a2, uint32_t *b1, uint32_t *b2, uint32_t *o) {
    while (b1 != b2) {
        // assert(0 <= *b1 && *b1 < a2 - a1)
        *o++ = a1[*b1++];
    }
}

uint32_t* foo_fx(uint32_t *a1, uint32_t *a2, uint32_t *b1, uint32_t *b2, const uint32_t b_offset, uint32_t *o) {
    while (b1 != b2) {
        // assert(0 <= *b1 && *b1 < a2 - a1)
        *o++ = a1[b_offset+(*b1++)];
    }
    return o;
}

uint32_t data[max_n], index[max_q], result[max_q];
std::pair<uint32_t, uint32_t[max_q / 8]>index_fx[16];

}

int main() {
    uint32_t seed = 0;
    auto rng = [&seed]() { return seed = seed * 9301 + 49297; };
    std::generate_n(data, max_n, rng);
    //std::generate_n(index, max_q, [rng]() { return rng() % max_n; });
    for (size_t i = 0; i < max_q;++i) {
        const uint32_t idx = rng() % max_n;
        const uint32_t bucket = idx >> 20;
        assert(bucket < 16);
        index_fx[bucket].second[index_fx[bucket].first] = idx % (1 << 20);
        index_fx[bucket].first++;
        assert((1 << 20)*bucket + index_fx[bucket].second[index_fx[bucket].first - 1] == idx);
    }
    auto t1 = std::chrono::high_resolution_clock::now();
    //foo(data, data + max_n, index, index + max_q, result);

    uint32_t* result_begin = result;
    for (int i = 0; i < 16; ++i) {
        result_begin = foo_fx(data, data + max_n, index_fx[i].second, index_fx[i].second + index_fx[i].first, (1<<20)*i, result_begin);
    }
    auto t2 = std::chrono::high_resolution_clock::now();
    std::cout << std::chrono::duration<double>(t2 - t1).count() << std::endl;
    std::cout << std::accumulate(result, result + max_q, 0ull) << std::endl;
}