c++ - 缓存未命中对矩阵乘法时间的影响

Question

我正在尝试进行矩阵乘法，从较小的矩阵大小开始并逐渐增加它，希望观察一旦矩阵不再适合缓存，时间会如何突然变化。但令我失望的是，在显然相同的函数之后，我总是得到一个非常平滑的图形。我尝试从小到 4x4 的矩阵大小开始，然后逐渐增加它直到 3400x3400 等于 11MB 整数值，但我仍然看不到时间函数的变化。可能是我在这里遗漏了一些关键点。任何帮助将不胜感激。

这是我的 C++ 代码:

long long clock_time() {
    struct timespec tp;
    clock_gettime(CLOCK_REALTIME, &tp);
    return (long long)(tp.tv_nsec + (long long)tp.tv_sec * 1000000000ll);
}

int main()
{
    for(int matrix_size = 100; matrix_size < 3500; matrix_size += 100)
    {
        int *A = new int[matrix_size*matrix_size];
        int *B = new int[matrix_size*matrix_size];
        int *C = new int[matrix_size*matrix_size];

        long long start = clock_time();

        for(int i = 0; i < matrix_size; ++i)
            for(int j = 0; j < matrix_size; ++j)
                for(int k = 0; k < matrix_size; ++k)
                {
                    C[i + j*matrix_size] = A[i + k*matrix_size] * B[k + j*matrix_size];
                }

        long long end = clock_time();
        long long totalTime = (end - start);

        std::cout << matrix_size << "," << totalTime << std::endl;

        delete[] A;
        delete[] B;
        delete[] C;
    }

    std::cout << "done" ;


    return 0;
}

这是我得到的数据的示例图:

详细数据可见https://docs.google.com/spreadsheets/d/1Xtri8w2sLZLQE0566Raducg7G2L4GLqNYIvP4nrp2t8/edit?usp=sharing

更新:根据 Zheyuan 和 Frank 的建议，我没有用值 i+j 初始化矩阵并将时间除以 2*N^3

for(int i = 0; i < matrix_size; i++)
{
    for(int j = 0; j < matrix_size; j++)
    {
        A[i + j * matrix_size] = i+j;
        B[i + j * matrix_size] = i+j;
        B[i + j * matrix_size] = i+j;
    }

}

结果如下:

更新 2:交换 i 和 j 循环后:

Answer 1

好吧，您一定会观察到计时的三次曲线。假设您使用两个 N * N 方阵，则矩阵乘法具有复杂性，或 2 * N ^ 3) 的浮点运算 (FLOP) 数量。随着 N 的增加，FLOP 的增加主导了时间的增加，您将不容易观察到延迟问题。

如果你想研究延迟的纯粹影响，你应该通过 FLOP 的数量“规范化”你的时间:

measured time / (2 * N ^ 3)

或者:

(2 * N ^ 3) / measured time

前者是每次 FLOP 花费的平均时间，而后者是每秒 FLOP，在文献中通常称为 FLOPs。 FLOPs 是性能的主要指标(至少对于科学计算而言)。预计随着 N 的增加，前一个指标会上升(延迟增加)，而后一个指标会下降(性能下降)。

很抱歉，我不会写 C++，所以无法修改您的代码(但这很容易，因为您只需要用 2 * N ^ 3 进行额外的除法)。我曾经用 C 代码做过同样的实验，这是我在 Intel Core 2 Duo 上的结果。请注意，我报告的是 MFLOP，即 10 ^ 6 FLOP。该图实际上是用 R 软件制作的。

---

我的上述观察实际上是假设您的其他一切都是正确的。但实际上，似乎并非如此。

首先，矩阵乘法是:

C[i + j*matrix_size] += A[i + k*matrix_size] * B[k + j*matrix_size];

注意 += 而不是 =。

其次，您的循环嵌套设计不当。您正在执行矩阵乘法 C = A * B，其中所有矩阵都存储在列主要顺序，因此您应该注意循环嵌套顺序以确保您始终有步幅-1 访问最内层循环。众所周知，j-k-i 循环嵌套在这种情况下是最优的。因此，请考虑以下事项:

for(int j = 0; j < matrix_size; ++j)
    for(int k = 0; k < matrix_size; ++k)
        for(int i = 0; i < matrix_size; ++i)
        {
            C[i + j*matrix_size] += A[i + k*matrix_size] * B[k + j*matrix_size];
        }

第三，您从矩阵大小 100 * 100 开始，它已经在 L1 缓存之外，主要是 64KB。我建议您从 N = 24 开始。一些文献表明，N = 60 是这种缓存的大致边界值。

第四，您需要重复多次乘法以消除测量偏差。目前，对于每次试验 N(或代码中的 matrix_size)，您执行一次乘法并测量时间。这是不准确的。对于小的 N 你会得到假的时间。如何重复 (1000/N + 1) ^ 3 次？

• 当N很小的时候，你重复了很多次；
• 随着 N 越来越接近 1000，您重复的次数会减少；
• 当 N > 1000 时，您实际上只执行一次乘法。

当然，不要忘记您需要将测量的时间除以重复的次数。

---

当然还有其他可以优化代码的地方，比如使用常量寄存器和消除地址计算中的整数乘法，但它们不太重要，因此没有涉及。数组的初始化也被跳过，因为它已经在 Frank 的回答中提出。