CPU缓存如何影响C程序的性能

Question

我试图更多地了解 CPU 缓存如何影响性能。作为一个简单的测试，我将矩阵第一列的值与不同数量的总列数相加。

// compiled with: gcc -Wall -Wextra -Ofast -march=native cache.c
// tested with: for n in {1..100}; do ./a.out $n; done | tee out.csv
#include <assert.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

double sum_column(uint64_t ni, uint64_t nj, double const data[ni][nj])
{
    double sum = 0.0;
    for (uint64_t i = 0; i < ni; ++i) {
        sum += data[i][0];
    }
    return sum;
}

int compare(void const* _a, void const* _b)
{
    double const a = *((double*)_a);
    double const b = *((double*)_b);
    return (a > b) - (a < b);
}

int main(int argc, char** argv)
{
    // set sizes
    assert(argc == 2);
    uint64_t const iter_max = 101;
    uint64_t const ni       = 1000000;
    uint64_t const nj       = strtol(argv[1], 0, 10);

    // initialize data
    double(*data)[nj] = calloc(ni, sizeof(*data));
    for (uint64_t i = 0; i < ni; ++i) {
        for (uint64_t j = 0; j < nj; ++j) {
            data[i][j] = rand() / (double)RAND_MAX;
        }
    }

    // test performance
    double* dt        = calloc(iter_max, sizeof(*dt));
    double const sum0 = sum_column(ni, nj, data);
    for (uint64_t iter = 0; iter < iter_max; ++iter) {
        clock_t const t_start = clock();
        double const sum      = sum_column(ni, nj, data);
        clock_t const t_stop  = clock();
        assert(sum == sum0);
        dt[iter] = (t_stop - t_start) / (double)CLOCKS_PER_SEC;
    }

    // sort dt
    qsort(dt, iter_max, sizeof(*dt), compare);

    // compute mean dt
    double dt_mean = 0.0;
    for (uint64_t iter = 0; iter < iter_max; ++iter) {
        dt_mean += dt[iter];
    }
    dt_mean /= iter_max;

    // print results
    printf("%2lu %.8e %.8e %.8e %.8e\n", nj, dt[iter_max / 2], dt_mean, dt[0],
        dt[iter_max - 1]);

    // free memory
    free(data);
}

然而，结果并不完全符合我的预期：

据我了解，当CPU从 加载值时data，它还会将以下一些值放入data缓存中。确切的数字取决于缓存行大小（我的机器上为 64 字节）。这可以解释为什么随着nj求解时间的增加，求解时间首先线性增加，并在某个值处趋于平稳。如果nj == 1，则一次加载会将接下来的 7 个值放入缓存中，因此我们只需每隔 8 个值从 RAM 加载一次。如果nj == 2，遵循相同的逻辑，我们需要每第 4 个值访问一次 RAM。在达到一定大小后，我们将必须访问每个值的 RAM，这应该会导致性能趋于平稳。我的猜测是，为什么该图的线性部分比 4 更进一步，因为实际上这里有多个级别的缓存在工作，并且值最终进入这些缓存的方式比我在这里解释的要复杂一些。

我无法解释的是为什么这些性能峰值出现在 16 的倍数处。

在思考了这个问题之后，我决定检查一下较高值是否也会发生这种情况nj：

事实上，确实如此。而且，还有更多：为什么性能在 ~250 之后再次提高？

有人可以向我解释一下，或者给我指出一些适当的参考资料，为什么会出现这些峰值，以及为什么性能会随着 的更高值而增加nj。

如果您想亲自尝试代码，为了方便起见，我还将附上我的绘图脚本：

import numpy as np
import matplotlib.pyplot as plt

data = np.genfromtxt("out.csv")
data[:,1:] /= data[0,1]

dy = np.diff(data[:,2]) / np.diff(data[:,0])
for i in range(len(dy) - 1):
    if dy[i] - dy[i + 1] > (dy.max() - dy.min()) / 2:
        plt.axvline(data[i + 1,0], color='gray', linestyle='--')
        plt.text(data[i + 1,0], 1.5 * data[0,3], f"{int(data[i + 1,0])}",
                 rotation=0, ha="center", va="center",
                 bbox=dict(boxstyle="round", ec='gray', fc='w'))

plt.fill_between(data[:,0], data[:,3], data[:,4], color='gray')
plt.plot(data[:,0], data[:,1], label="median")
plt.plot(data[:,0], data[:,2], label="mean")
plt.legend(loc="upper left")
plt.xlabel("nj")
plt.ylabel("dt / dt$_0$")
plt.savefig("out.pdf")

Answer 1

这些图显示了几种复杂的低级效应的组合（主要是缓存垃圾和预取问题）。我假设目标平台是一款主流现代处理器，具有 64 字节缓存行（通常是 x86）。

我可以在我的 i5-9600KF 处理器上重现该问题。这是结果图：

首先，当nj较小时，获取的地址之间的间隙（即步幅）较小，并且缓存线的使用相对有效。例如，当 时nj = 1，访问是连续的。在这种情况下，处理器可以有效地从 DRAM预取缓存线，从而隐藏其高延迟。由于许多连续项共享相同的缓存行，因此还具有良好的空间缓存局部性。当 时nj=2，仅使用缓存行的一半值。这意味着对于相同数量的操作，请求的缓存行数量是两倍。话虽这么说，由于将两个浮点数相加导致计算密集型代码的延迟相对较高，所以时间并没有长很多。您可以展开循环4 次并使用 4 个不同的总和变量，以便（主流现代）处理器可以并行添加多个值。请注意，大多数处理器还可以在每个周期从缓存加载多个值。每 2 个周期请求nj = 4一个新的缓存行（因为 adouble占用 8 个字节）。结果，内存吞吐量可能变得如此之大，以至于计算变得受内存限制。人们可能期望时间是稳定的，nj >= 8因为请求的缓存行的数量应该相同，但实际上，处理器预取多个连续的缓存行，这样就不用支付在这种情况下巨大的 DRAM 延迟的开销。预取缓存行的数量通常在 2 到 4 之间（据我所知，当步幅大于 512 时，这种预取策略在 Intel 处理器上被禁用，因此当nj >= 64。这解释了为什么时序急剧增加，nj < 32并且它们变得相对稳定，但有32 <= nj <= 256例外）对于峰值。

nj当是 16 的倍数时出现的常规峰值是由于称为缓存抖动的复杂缓存效应造成的。现代缓存是N 路关联的，N 通常在 4 到 16 之间。例如，以下是我的 i5-9600KF 处理器的统计数据：

Cache 0: L1 data cache,        line size 64,  8-ways,    64 sets, size 32k 
Cache 1: L1 instruction cache, line size 64,  8-ways,    64 sets, size 32k 
Cache 2: L2 unified cache,     line size 64,  4-ways,  1024 sets, size 256k 
Cache 3: L3 unified cache,     line size 64, 12-ways, 12288 sets, size 9216k 

这意味着，如果 的话，从 DRAM 中分别使用地址 A1 和 A2 获取的两个值可能会导致我的 L1 缓存中发生冲突(A1 % 32768) / 64 == (A2 % 32768) / 64。在这种情况下，处理器需要从一组缓存行中选择要替换N=8的缓存行。缓存替换策略有很多，但没有一个是完美的。因此，某些有用的高速缓存行有时会被过早逐出，从而导致稍后需要额外的高速缓存未命中。在病理情况下，许多 DRAM 位置可能会竞争相同的高速缓存行，从而导致过多的高速缓存未命中。有关此内容的更多信息也可以在这篇文章中找到。

关于nj步幅，L1缓存中可以有效使用的缓存行的数量是有限的。例如，如果所有获取的值具有相同的地址模数和高速缓存大小，则实际上只能使用 N 个高速缓存行（即，对于我的处理器而言为 8 个）来存储所有值。可用的缓存行较少是一个大问题，因为预取器需要缓存中相当大的空间，以便存储稍后需要的许多缓存行。并发读取数量越少，内存吞吐量越低。此处尤其如此，因为从 DRAM 获取 1 个缓存行的延迟约为几十纳秒（例如 ~70 ns），而其带宽约为数十 GiB/s（例如 ~40 GiB/s）：应同时获取高速缓存行（例如~40），以隐藏延迟并使 DRAM 饱和。

以下是关于 的值的 L1 缓存中实际可以使用的缓存行数的模拟nj：

 nj  #cache-lines
  1      512
  2      512
  3      512
  4      512
  5      512
  6      512
  7      512
  8      512
  9      512
 10      512
 11      512
 12      512
 13      512
 14      512
 15      512
 16      256    <----
 17      512
 18      512
 19      512
 20      512
 21      512
 22      512
 23      512
 24      512
 25      512
 26      512
 27      512
 28      512
 29      512
 30      512
 31      512
 32      128    <----
 33      512
 34      512
 35      512
 36      512
 37      512
 38      512
 39      512
 40      512
 41      512
 42      512
 43      512
 44      512
 45      512
 46      512
 47      512
 48      256    <----
 49      512
 50      512
 51      512
 52      512
 53      512
 54      512
 55      512
 56      512
 57      512
 58      512
 59      512
 60      512
 61      512
 62      512
 63      512
 64       64    <----
==============
 80      256
 96      128
112      256
128       32
144      256
160      128
176      256
192       64
208      256
224      128
240      256
256       16
384       32
512        8
1024       4

我们可以看到，当 是 16 的倍数时，可用缓存行的数量较小。nj在这种情况下，预加载器会将数据预加载到可能被后续提取（同时完成）提前逐出的缓存行中。当可用高速缓存行的数量较少时，代码中执行的加载指令更有可能导致高速缓存未命中。当发生缓存未命中时，需要再次从 L2 甚至 L3 获取该值，从而导致执行速度变慢。请注意，L2 缓存也受到相同的影响，但由于其较大而不太明显。现代 x86 处理器的 L3 缓存利用散列来更好地分配事物，以减少固定步幅的冲突（至少在 Intel 处理器上，当然在 AMD 上也是如此，尽管 AFAIK 这没有记录在案）。

以下是我的机器上一些峰值的计时：

  32 4.63600000e-03 4.62298020e-03 4.06400000e-03 4.97300000e-03
  48 4.95800000e-03 4.96994059e-03 4.60400000e-03 5.59800000e-03
  64 5.01600000e-03 5.00479208e-03 4.26900000e-03 5.33100000e-03
  96 4.99300000e-03 5.02284158e-03 4.94700000e-03 5.29700000e-03
 128 5.23300000e-03 5.26405941e-03 4.93200000e-03 5.85100000e-03
 192 4.76900000e-03 4.78833663e-03 4.60100000e-03 5.01600000e-03
 256 5.78500000e-03 5.81666337e-03 5.77600000e-03 6.35300000e-03
 384 5.25900000e-03 5.32504950e-03 5.22800000e-03 6.75800000e-03
 512 5.02700000e-03 5.05165347e-03 5.02100000e-03 5.34400000e-03
1024 5.29200000e-03 5.33059406e-03 5.28700000e-03 5.65700000e-03

正如预期的那样，在可用缓存行数量少得多的情况下，实践中的时序总体上会更大。然而，当 时nj >= 512，结果令人惊讶，因为它们比其他方法快得多。在这种情况下，可用高速缓存行的数量等于关联方式的数量 (N)。我的猜测是，这是因为英特尔处理器肯定会检测到这种病态情况并优化预取，以减少缓存未命中的数量（使用行填充缓冲区绕过 L1 缓存 - 见下文）。

最后，对于大nj步幅，较大的nj应该会导致较高的开销，这主要是由于转换后备缓冲区（TLB）：有更多的页面地址需要较大的转换nj，并且 TLB 条目的数量是有限的。事实上，这是我在我的机器上可以观察到的情况：与目标平台不同，时间往往以非常稳定的方式缓慢增加。

我还无法真正解释这种非常奇怪的行为。以下是一些疯狂的猜测：

• 当 较大时，操作系统可能倾向于使用更多的大页nj（以便减少 TLB 的开销），因为分配了更宽的块。这可能会导致预取器有更多的并发性，因为据我所知它不能跨越页面边界。您可以尝试检查分配的（透明）大页面的数量（通过在 Linux 中查看AnonHugePages）/proc/meminfo或在这种情况下强制使用它们（使用显式memmap），或者可能通过禁用它们。我的系统似乎使用 2 MiB 透明大页面，与nj值无关。
• 如果目标体系结构是 NUMA 体系结构（例如，新的 AMD 处理器或具有多个处理器且拥有自己的内存的服务器），则操作系统可以分配物理存储在另一个 NUMA 节点上的页面，因为当前 NUMA 节点上的可用空间较少。由于吞吐量更大，这可能会带来更高的性能（尽管延迟更高）。您可以在 Linux 上控制此策略numactl，以便强制进行本地分配。

有关此主题的更多信息，请阅读精彩文档《每个程序员都应该了解内存》。此外，这里还提供了一篇关于 x86 缓存在实践中如何工作的非常好的文章。

去除峰

要消除 x86 处理器上由于缓存垃圾造成的峰值，您可以使用非临时软件预取指令，以便可以在非临时缓存结构中提取缓存行并将其提取到靠近处理器的位置，这样不会导致处理器中的缓存垃圾。 L1（如果可能）。这种缓存结构通常是Intel 处理器上的行填充缓冲区(LFB) 和 AMD Zen 处理器上的（等效）未命中地址缓冲区 (MAB)。有关非临时指令和 LFB 的更多信息，请阅读这篇文章和这篇文章。nj下面是修改后的代码，其中还包括循环展开优化，以在较小时加快代码速度：

double sum_column(uint64_t ni, uint64_t nj, double* const data)
{
    double sum0 = 0.0;
    double sum1 = 0.0;
    double sum2 = 0.0;
    double sum3 = 0.0;

    if(nj % 16 == 0)
    {
        // Cache-bypassing prefetch to avoid cache trashing
        const size_t distance = 12;
        for (uint64_t i = 0; i < ni; ++i) {
            _mm_prefetch(&data[(i+distance)*nj+0], _MM_HINT_NTA);
            sum0 += data[i*nj+0];
        }
    }
    else
    {
        // Unrolling is much better for small strides
        for (uint64_t i = 0; i < ni; i+=4) {
            sum0 += data[(i+0)*nj+0];
            sum1 += data[(i+1)*nj+0];
            sum2 += data[(i+2)*nj+0];
            sum3 += data[(i+3)*nj+0];
        }
    }
    
    return sum0 + sum1 + sum2 + sum3;
}

这是修改后的代码的结果：

我们可以看到峰值不再出现在时间中。我们还可以看到，由于值小了dt0大约 4 倍（由于循环展开），所以值要大得多。

请注意，在实践中，这种方法并不能避免二级缓存中的缓存垃圾（至少在英特尔处理器上）。这意味着在我的机器上，效果仍然存在，nj步幅为 512 (4 KiB) 的倍数（实际上比以前慢，尤其是在 时nj >= 2048）。在 x86 处理器上停止预取可能是个好主意(nj%512) == 0 && nj >= 512。效果 AFAIK，没有办法解决这个问题。话虽这么说，对非常大的数据结构执行如此大的跨步访问是一个非常糟糕的主意。

请注意，distance应仔细选择，因为早期预取可能会导致缓存行在实际使用之前被逐出（因此需要再次获取），而后期预取则没有多大用处。我认为使用接近 LFB/MAB 中条目数量的值是一个好主意（例如，Skylake/KabyLake/CannonLake 上为 12，Zen-2 上为 22）。