用于缓存敏感矩阵乘法的多级平铺答案

【问题标题】：Multilevel tiling for cache sensitive matrix multiplication用于缓存敏感矩阵乘法的多级平铺
【发布时间】：2021-05-06 15:04:55
【问题描述】：

为了好玩，我正在编写自己的 GeMM 子程序。我已经设法使用 AVX256 内核在 L1 缓存上实现了平铺版本。我知道有一些循环不变量我可以从下面提升。

template<size_t N, size_t P, size_t M>
void intrinsic_multiply(double* A, double* B, void* C) {
    alignas(32) double _A[PACK_SIZE * PACK_SIZE];
    alignas(32) double _B[PACK_SIZE * PACK_SIZE];
    double* _C = (double*)(C);

    constexpr size_t N_rem = N%PACK_SIZE;
    constexpr size_t M_rem = M%PACK_SIZE;
    constexpr size_t P_rem = P%PACK_SIZE;

    constexpr size_t N_spill = N-N_rem+(PACK_SIZE)*(size_t)(N_rem!=0);
    constexpr size_t M_spill = M-M_rem+(PACK_SIZE)*(size_t)(M_rem!=0);
    constexpr size_t P_spill = P-P_rem+(PACK_SIZE)*(size_t)(P_rem!=0);

    for(size_t i = 0; i < N_spill; i += PACK_SIZE) {
        for(size_t j = 0; j < M_spill; j += PACK_SIZE) {
            for(size_t k = 0; k < P_spill; k += PACK_SIZE) {
                pad_cols<N>(A + i + k * N, _A, ((size_t)(i+PACK_SIZE>N))*N_rem, ((size_t)(k+PACK_SIZE>P))*P_rem);
                pad_cols<P>(B + k + j * P, _B, ((size_t)(k+PACK_SIZE>P))*P_rem, ((size_t)(j+PACK_SIZE>M))*M_rem);
                macro_kernal_intrinsic<N_spill>(_A, _B, _C + i + (j * N_spill));
            }
        }
    }
}

我很难实现多级缓存平铺，因为缓存不是彼此的倍数。每个缓存的步长估计值计算如下，CACHE_SIZEs 以字节为单位。

static constexpr size_t L3_CACHESIZE = 6291456;
static size_t L3_STRIDE_SIZE = (size_t)floor((sqrt(L3_CACHESIZE/sizeof(double))));

static constexpr size_t L2_CACHESIZE = 262144;
static size_t L2_STRIDE_SIZE = (size_t)floor((sqrt(L2_CACHESIZE/sizeof(double))));

static constexpr size_t L1_CACHESIZE = 32768;
static size_t L1_STRIDE_SIZE = (size_t)floor((sqrt(L1_CACHESIZE/sizeof(double))));

static constexpr size_t PACK_SIZE = 64;

这给出了 L1 步幅大小为 64 和 L2 步幅大小为 181。显然，它们不是彼此的倍数。我有两个选择 -

在每个 L2 迭代中仅适合 4 个 L1 块，从 0 到 63 到 127。这似乎是我正在利用我的 L2 缓存。
使用整个二级缓存并在最后一次迭代中用 0 填充 (64*3-180) 个元素。这会引入大量冗余操作，但只会将 L2 步长减小 1。
为下一次迭代预取一个小数块。
有一种我不知道的规范方法。

在实践中处理不是彼此倍数的块大小的最佳方法是什么？

编辑 - 响应 MSalters：

我跑了

nm -AP -t d --demangle ./src/example-GeMM | grep "intrinsic"

给了

./src/example-GeMM: void intrinsic_multiply<2048ul, 3136ul, 2240ul>(double*, double*, void*) W 17081 434
./src/example-GeMM: void macro_kernal_intrinsic<2048ul>(double*, double*, double*) W 20409 234
./src/example-GeMM: void micro_kernal_intrinsic<2048ul>(double*, double*, double*) W 21343 454

意味着相关代码段占用 (234+454+434)/262144

【问题讨论】：

_A 是保留名称，即未定义的行为。
您将无法使用“经典”（行或列专业）排序超过特定限制。所有“更快”的矩阵乘法算法都会将矩阵重新排序为一些更“方便”的排序。

标签： c++ matrix-multiplication cpu-cache tiling

【解决方案1】：

我的假设是 L1 被拆分，然后 32Kb 只是 L1d 缓存。 L2 将是统一的，这意味着您不应仅将其全部用于您的数据。

另外，我还不太担心循环提升。有了所有的constexpr，优化器就有了发现提升机的机会。一个问题是C 似乎是外变量。因为那是void*，它可以别名为size_t。如果你写了double* C，这不能别名size_t。这是一个示例，说明类型安全不仅有利于安全，而且有利于速度。

【讨论】：

我会更改_A。是的，L1 是从 64kb 缓存中分离出来的。感谢您提供有关别名的提示。 L3是共享缓存，每个处理器都有自己的L1和L2缓存。希望有人仍然可以就处理不同的缓存大小提出一些明确的建议。
@IvorDenham-Dyson：当我说 L2 是统一的时，我想你错过了我的意思。这并不意味着它在核心之间共享。这意味着它在代码和数据之间共享。因此，您“未充分利用”L2 的想法很好，因为它在 L2 中为代码留出了空间。
好点，不知道L2大体上是统一的。请查看我的编辑。

【解决方案2】：

请参阅以下论文（scholar.google.com 上的类似论文）。我认为这个很适合你的问题。（顺便说一句，我想你已经使用了这个“_mm256_mul_pd”。）无论如何，你会在互联网上找到pdf（不想在这里链接）。

GPU 上高效 GEMM 的协调平铺和批处理框架， X Li, Y Liang, S Yan, L Jia, Y Li - 24th Symposium on …, 2019 - dl.acm.org

“最根本的问题是平铺、批处理以及它们的协同交互。平铺意味着将每个GEMM平铺成许多平铺。我们允许不同的GEMM有不同的平铺策略，而不是共享一个统一的平铺策略。如何统一不同的平铺策略进入单个内核是一项挑战。”

"给定一个大小为 M×N×K 的 GEMM，将 C 矩阵划分为多个大小为 BY×BX 的图块。C 的每个图块需要访问大小为 BY×K 的 A 矩阵的一整行部分和一个图 1(a) 中大小为 K×BX 的 B 矩阵的整列段，但是 A 的整行带和 B 的整列带太大，无法容纳在共享存储器和寄存器文件中。内存，工作量沿K维划分为如图1（b）所示的许多段。A的行部分的每个部分称为A tile，大小为BY×BK，列部分的每个部分B的大小称为B瓦片，大小为BK×BX。最终的结果可以通过将每个段沿K维的部分结果累加得到。"

评论：我喜欢考虑不同解决方案的想法 - 或者让它“优化”以找到最佳解决方案。（比如：你为什么要找出哪个瓷砖尺寸表现最好，让计算机找出来。嗯，很明显，编写一个瓷砖尺寸的可变设置更费力）。好吧，这是研究（我只知道您最好以常规顺序存储和访问 RAM 中的数据 - 否则您会看到“缓存未命中”，这可能会导致巨大的性能损失。另一个想法：其他一些进程可能会喜欢也住在缓存中。我不确定你是如何考虑到这一点的——没有提到。）

【讨论】：