具有独特矩阵转置问题的二维分块
2D blocking with unique matrix transpose problem
我有以 3 阶张量形式组织的 struct complex {double real = 0.0; double imag = 0.0;}; 类型的复数值数据。底层容器具有与内存页面边界对齐的连续内存布局。
张量的自然 'slicing' 方向是沿方向 1。这意味着缓存行按方向 3、2 和最后 1 的顺序延伸。换句话说，索引函数看起来像这样：(i, j, k) -> i * N2 * N3 + j * N3 + k。
我需要沿方向 2 转置切片。在上面的第一张图像中，红色矩形是我希望转置的张量切片。
我的 C++ 代码如下所示：
for (auto vslice_counter = std::size_t{}; vslice_counter < n2; ++vslice_counter)
{
// blocked loop
for (auto bi1 = std::size_t{}; bi1 < n1; bi1 += block_size)
{
for (auto bi3 = std::size_t{}; bi3 < n3; bi3 += block_size)
{
for (auto i1 = std::size_t{}; i1 < block_size; ++i1)
{
for (auto i3 = std::size_t{}; i3 < block_size; ++i3)
{
const auto i1_block = bi1 + i1;
const auto i3_block = bi3 + i3;
tens_tr(i3_block, vslice_counter, i1_block) = tens(i1_block, vslice_counter, i3_block);
}
}
}
}
}
用于测试的机器：双路 Intel(R) Xeon(R) Platinum 8168 with
- 24 核 @ 2.70 GHz 和
- L1、L2 和 L3 缓存大小分别为 32 kB、1 MB 和 33 MB。
我绘制了此函数 w.r.t 块大小的性能图，但令我惊讶的是没有变化！事实上，朴素的实现与这个实现一样好。
问题：在这种情况下，2D 分块无助于提高性能是否有原因？
编辑：
附加信息：
- 使用的编译器是Intel C++编译器。
block_size是输入参数，不是编译时常量。- 在测试期间，我将
block_size 的值从 4 更改为 256。这分别转换为 256 B 和 1 MB 的块，因为这些块是边长 block_size 的正方形。我的目标是用块填充 L1 and/or L2 缓存。
- 用于优化的编译标志：
-O3;-ffast-math;-march=native;-qopt-zmm-usage=high
TL;DR：2D 阻塞代码并不快，主要是因为 缓存垃圾处理 （还有一点是由于预取）。简单的实现和 2D 阻塞都是低效的。 中描述了类似的效果。您可以使用小型 临时数组 来缓解此问题，并使用 register-tiling 优化来提高性能。额外的填充也可能有帮助。
CPU 缓存
首先，现代 CPU 缓存是 set-associative cache。 m-way 关联缓存可以看作是一个 n x m 矩阵，其中 n 集包含一块 m 缓存行。内存块首先映射到一个集合上，然后放入目标集合的任何 m 缓存行（关于缓存替换策略）。当算法（如转置）执行跨步内存访问时，处理器经常访问同一组，并且可以使用的目标缓存行数要少得多。对于 2 的大幂次（或更准确地说可被 2 的大幂次整除）的步幅尤其如此。可以使用模块化算法或基本模拟来计算可以访问的可能缓存行的数量。在最坏的情况下，所有访问都在同一缓存行集（包含 m 缓存行）上完成。在这种情况下，每次访问都会导致处理器使用通常不完美的替换策略逐出集合中的一个缓存行以在其中加载一个新缓存行。您的 Intel(R) Xeon(R) Platinum 8168 处理器具有以下缓存属性：
Level Size (KiB/core) Associativity # Sets
-------------------------------------------------------
L1D Cache 32 8-way 64
L2 Cache 1024 16-way 1024
L3 Cache 1408 11-way 2048
* All cache have 64-byte cache lines
这意味着您执行的访问步幅可被 4096 字节整除（即 256 个复数），然后所有访问都将映射到 L1D 缓存中的相同 L1D 缓存集，从而导致 冲突未命中 因为只能同时加载 8 个缓存行。这会导致称为 缓存垃圾处理 的效果大大降低性能。 post 中有更完整的解释：.
缓存对换位的影响
您提到 block_size 可以达到 256，因此 n1 和 n3 可以被 256 整除，因为提供的代码已经隐含假设 n1 和n3 可以被 block_size 整除，我希望 n2 也可以被 256 整除。这意味着沿维度 3 的线的大小为 p * 256 * (2 * 8) = 4096p bytes = 64p cache lines，其中 p = n3 / block_size。这意味着所有行的第 i 个项目将映射到完全相同的 L1D 缓存集。
由于在维度 1 和维度 3 之间进行了换位，这意味着行与行之间的 space 更大。两个后续行的第 i 个项目之间的差距是 G = 64 * p * n2 个缓存行。假设 n2 可以被 16 整除，那么 G = 1024 * p * q 缓存行 q = n2 / 16.
有这样的差距是个大问题。实际上，您的算法 reads/writes 到同一块中多行的许多第 i 项。因此，此类访问将映射到 L1D 和 L2 缓存中的同一组，从而导致缓存垃圾。如果 block_size > 16，缓存行将几乎系统地重新加载到 L2（4 次）。我预计 L3 在这种情况下不会有太大帮助，因为它主要是为内核之间的共享数据而设计的，它的延迟非常大（50-70 个周期）并且 p * q 肯定可以被 2 的幂整除。由于 缺乏并发性 （即可以同时预取的可用缓存行），处理器无法缓解延迟。这导致带宽被浪费，更不用说 non-contiguous 访问已经降低了吞吐量。这种效果应该已经用两个 block_size 值的较小幂可见，如先前相关 post（上面链接）所示。
预取的影响
英特尔 Skylake SP 处理器，如您的 prefetch 在这种情况下每次访问同时至少 2 个缓存行（128 字节）。这意味着 block_size < 8 不足以完全使用 tens 的预取缓存行。结果，a 太小 block_size 由于预取而浪费带宽，太大也浪费带宽，但由于缓存垃圾 。我希望最好的 block_size 接近 8。
如何解决这个问题
一个解决方案是将每个块存储在一个小的临时二维数组中，然后转置它，然后写入它。乍一看，好像是因为访问内存多了会变慢，但一般来说速度会快很多。事实上，如果块大小相当小（例如 <=32），那么小的临时数组可以完全适合 L1 缓存，因此在转置过程中不会受到任何缓存垃圾的影响。该块可以同样有效地读取，但可以更有效地存储它（即更连续的访问）。
添加另一个阻塞级别应该有助于提高性能，因为 L2 缓存可以更有效地使用（例如 block_size 设置 128~256）。 Lebesgue curve可以用o 实施快速 缓存遗忘 算法，尽管它会使代码更复杂。
另一个优化包括添加另一个阻塞级别来执行称为 register-tiling 的优化。这个想法是使用 2 个嵌套循环操作一个带有小 compile-time 常量 的图块，以便编译器 展开循环 并生成更好的说明。例如，对于大小为 4x4 的图块，这使编译器能够生成以下代码（参见 Godbolt）：
..B3.7:
vmovupd xmm0, XMMWORD PTR [rdx+r8]
vmovupd XMMWORD PTR [r15+rdi], xmm0
inc r14
vmovupd xmm1, XMMWORD PTR [16+rdx+r8]
vmovupd XMMWORD PTR [r15+r10], xmm1
vmovupd xmm2, XMMWORD PTR [32+rdx+r8]
vmovupd XMMWORD PTR [r15+r12], xmm2
vmovupd xmm3, XMMWORD PTR [48+rdx+r8]
vmovupd XMMWORD PTR [r15+r13], xmm3
vmovupd xmm4, XMMWORD PTR [rdx+r9]
vmovupd XMMWORD PTR [16+r15+rdi], xmm4
vmovupd xmm5, XMMWORD PTR [16+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r10], xmm5
vmovupd xmm6, XMMWORD PTR [32+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r12], xmm6
vmovupd xmm7, XMMWORD PTR [48+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r13], xmm7
vmovupd xmm8, XMMWORD PTR [rdx+r11]
vmovupd XMMWORD PTR [32+r15+rdi], xmm8
vmovupd xmm9, XMMWORD PTR [16+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r10], xmm9
vmovupd xmm10, XMMWORD PTR [32+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r12], xmm10
vmovupd xmm11, XMMWORD PTR [48+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r13], xmm11
vmovupd xmm12, XMMWORD PTR [rdx+rbp]
vmovupd XMMWORD PTR [48+r15+rdi], xmm12
vmovupd xmm13, XMMWORD PTR [16+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r10], xmm13
vmovupd xmm14, XMMWORD PTR [32+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r12], xmm14
vmovupd xmm15, XMMWORD PTR [48+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r13], xmm15
add r15, rsi
add rdx, 64
cmp r14, rbx
jb ..B3.7
代替这个（重复8次以上）：
..B2.12:
vmovupd xmm0, XMMWORD PTR [rsi+r14]
vmovupd XMMWORD PTR [rbx+r15], xmm0
inc rax
vmovupd xmm1, XMMWORD PTR [16+rsi+r14]
vmovupd XMMWORD PTR [rbx+r13], xmm1
add rbx, r9
add rsi, 32
cmp rax, rcx
jb ..B2.12
最后，可以使用 AVX/AVX-2/AVX-512 内在函数为仅 x86-64 处理器实现更快的图块转置。
请注意，在每行的末尾添加一些填充，使它们不能被 的幂整除，这也应该显着有助于减少缓存垃圾。也就是说，一旦应用了上述优化，这可能就没有用了。
我有以 3 阶张量形式组织的 struct complex {double real = 0.0; double imag = 0.0;}; 类型的复数值数据。底层容器具有与内存页面边界对齐的连续内存布局。
张量的自然 'slicing' 方向是沿方向 1。这意味着缓存行按方向 3、2 和最后 1 的顺序延伸。换句话说，索引函数看起来像这样：(i, j, k) -> i * N2 * N3 + j * N3 + k。
我需要沿方向 2 转置切片。在上面的第一张图像中，红色矩形是我希望转置的张量切片。
我的 C++ 代码如下所示：
for (auto vslice_counter = std::size_t{}; vslice_counter < n2; ++vslice_counter)
{
// blocked loop
for (auto bi1 = std::size_t{}; bi1 < n1; bi1 += block_size)
{
for (auto bi3 = std::size_t{}; bi3 < n3; bi3 += block_size)
{
for (auto i1 = std::size_t{}; i1 < block_size; ++i1)
{
for (auto i3 = std::size_t{}; i3 < block_size; ++i3)
{
const auto i1_block = bi1 + i1;
const auto i3_block = bi3 + i3;
tens_tr(i3_block, vslice_counter, i1_block) = tens(i1_block, vslice_counter, i3_block);
}
}
}
}
}
用于测试的机器：双路 Intel(R) Xeon(R) Platinum 8168 with
- 24 核 @ 2.70 GHz 和
- L1、L2 和 L3 缓存大小分别为 32 kB、1 MB 和 33 MB。
我绘制了此函数 w.r.t 块大小的性能图，但令我惊讶的是没有变化！事实上，朴素的实现与这个实现一样好。
问题：在这种情况下，2D 分块无助于提高性能是否有原因？
编辑：
附加信息：
- 使用的编译器是Intel C++编译器。
block_size是输入参数，不是编译时常量。- 在测试期间，我将
block_size的值从 4 更改为 256。这分别转换为 256 B 和 1 MB 的块，因为这些块是边长block_size的正方形。我的目标是用块填充 L1 and/or L2 缓存。 - 用于优化的编译标志：
-O3;-ffast-math;-march=native;-qopt-zmm-usage=high
TL;DR：2D 阻塞代码并不快，主要是因为 缓存垃圾处理 （还有一点是由于预取）。简单的实现和 2D 阻塞都是低效的。
CPU 缓存
首先，现代 CPU 缓存是 set-associative cache。 m-way 关联缓存可以看作是一个 n x m 矩阵，其中 n 集包含一块 m 缓存行。内存块首先映射到一个集合上，然后放入目标集合的任何 m 缓存行（关于缓存替换策略）。当算法（如转置）执行跨步内存访问时，处理器经常访问同一组，并且可以使用的目标缓存行数要少得多。对于 2 的大幂次（或更准确地说可被 2 的大幂次整除）的步幅尤其如此。可以使用模块化算法或基本模拟来计算可以访问的可能缓存行的数量。在最坏的情况下，所有访问都在同一缓存行集（包含 m 缓存行）上完成。在这种情况下，每次访问都会导致处理器使用通常不完美的替换策略逐出集合中的一个缓存行以在其中加载一个新缓存行。您的 Intel(R) Xeon(R) Platinum 8168 处理器具有以下缓存属性：
Level Size (KiB/core) Associativity # Sets
-------------------------------------------------------
L1D Cache 32 8-way 64
L2 Cache 1024 16-way 1024
L3 Cache 1408 11-way 2048
* All cache have 64-byte cache lines
这意味着您执行的访问步幅可被 4096 字节整除（即 256 个复数），然后所有访问都将映射到 L1D 缓存中的相同 L1D 缓存集，从而导致 冲突未命中 因为只能同时加载 8 个缓存行。这会导致称为 缓存垃圾处理 的效果大大降低性能。 post 中有更完整的解释：
缓存对换位的影响
您提到 block_size 可以达到 256，因此 n1 和 n3 可以被 256 整除，因为提供的代码已经隐含假设 n1 和n3 可以被 block_size 整除，我希望 n2 也可以被 256 整除。这意味着沿维度 3 的线的大小为 p * 256 * (2 * 8) = 4096p bytes = 64p cache lines，其中 p = n3 / block_size。这意味着所有行的第 i 个项目将映射到完全相同的 L1D 缓存集。
由于在维度 1 和维度 3 之间进行了换位，这意味着行与行之间的 space 更大。两个后续行的第 i 个项目之间的差距是 G = 64 * p * n2 个缓存行。假设 n2 可以被 16 整除，那么 G = 1024 * p * q 缓存行 q = n2 / 16.
有这样的差距是个大问题。实际上，您的算法 reads/writes 到同一块中多行的许多第 i 项。因此，此类访问将映射到 L1D 和 L2 缓存中的同一组，从而导致缓存垃圾。如果 block_size > 16，缓存行将几乎系统地重新加载到 L2（4 次）。我预计 L3 在这种情况下不会有太大帮助，因为它主要是为内核之间的共享数据而设计的，它的延迟非常大（50-70 个周期）并且 p * q 肯定可以被 2 的幂整除。由于 缺乏并发性 （即可以同时预取的可用缓存行），处理器无法缓解延迟。这导致带宽被浪费，更不用说 non-contiguous 访问已经降低了吞吐量。这种效果应该已经用两个 block_size 值的较小幂可见，如先前相关 post（上面链接）所示。
预取的影响
英特尔 Skylake SP 处理器，如您的 prefetch 在这种情况下每次访问同时至少 2 个缓存行（128 字节）。这意味着 block_size < 8 不足以完全使用 tens 的预取缓存行。结果，a 太小 block_size 由于预取而浪费带宽，太大也浪费带宽，但由于缓存垃圾 。我希望最好的 block_size 接近 8。
如何解决这个问题
一个解决方案是将每个块存储在一个小的临时二维数组中，然后转置它，然后写入它。乍一看，好像是因为访问内存多了会变慢，但一般来说速度会快很多。事实上，如果块大小相当小（例如 <=32），那么小的临时数组可以完全适合 L1 缓存，因此在转置过程中不会受到任何缓存垃圾的影响。该块可以同样有效地读取，但可以更有效地存储它（即更连续的访问）。
添加另一个阻塞级别应该有助于提高性能，因为 L2 缓存可以更有效地使用（例如 block_size 设置 128~256）。 Lebesgue curve可以用o 实施快速 缓存遗忘 算法，尽管它会使代码更复杂。
另一个优化包括添加另一个阻塞级别来执行称为 register-tiling 的优化。这个想法是使用 2 个嵌套循环操作一个带有小 compile-time 常量 的图块，以便编译器 展开循环 并生成更好的说明。例如，对于大小为 4x4 的图块，这使编译器能够生成以下代码（参见 Godbolt）：
..B3.7:
vmovupd xmm0, XMMWORD PTR [rdx+r8]
vmovupd XMMWORD PTR [r15+rdi], xmm0
inc r14
vmovupd xmm1, XMMWORD PTR [16+rdx+r8]
vmovupd XMMWORD PTR [r15+r10], xmm1
vmovupd xmm2, XMMWORD PTR [32+rdx+r8]
vmovupd XMMWORD PTR [r15+r12], xmm2
vmovupd xmm3, XMMWORD PTR [48+rdx+r8]
vmovupd XMMWORD PTR [r15+r13], xmm3
vmovupd xmm4, XMMWORD PTR [rdx+r9]
vmovupd XMMWORD PTR [16+r15+rdi], xmm4
vmovupd xmm5, XMMWORD PTR [16+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r10], xmm5
vmovupd xmm6, XMMWORD PTR [32+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r12], xmm6
vmovupd xmm7, XMMWORD PTR [48+rdx+r9]
vmovupd XMMWORD PTR [16+r15+r13], xmm7
vmovupd xmm8, XMMWORD PTR [rdx+r11]
vmovupd XMMWORD PTR [32+r15+rdi], xmm8
vmovupd xmm9, XMMWORD PTR [16+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r10], xmm9
vmovupd xmm10, XMMWORD PTR [32+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r12], xmm10
vmovupd xmm11, XMMWORD PTR [48+rdx+r11]
vmovupd XMMWORD PTR [32+r15+r13], xmm11
vmovupd xmm12, XMMWORD PTR [rdx+rbp]
vmovupd XMMWORD PTR [48+r15+rdi], xmm12
vmovupd xmm13, XMMWORD PTR [16+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r10], xmm13
vmovupd xmm14, XMMWORD PTR [32+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r12], xmm14
vmovupd xmm15, XMMWORD PTR [48+rdx+rbp]
vmovupd XMMWORD PTR [48+r15+r13], xmm15
add r15, rsi
add rdx, 64
cmp r14, rbx
jb ..B3.7
代替这个（重复8次以上）：
..B2.12:
vmovupd xmm0, XMMWORD PTR [rsi+r14]
vmovupd XMMWORD PTR [rbx+r15], xmm0
inc rax
vmovupd xmm1, XMMWORD PTR [16+rsi+r14]
vmovupd XMMWORD PTR [rbx+r13], xmm1
add rbx, r9
add rsi, 32
cmp rax, rcx
jb ..B2.12
最后，可以使用 AVX/AVX-2/AVX-512 内在函数为仅 x86-64 处理器实现更快的图块转置。
请注意，在每行的末尾添加一些填充，使它们不能被 的幂整除，这也应该显着有助于减少缓存垃圾。也就是说，一旦应用了上述优化，这可能就没有用了。