为什么转置 512x512 的矩阵比转置 513x513 的矩阵慢得多？

解释来自 Agner Fog 在 C++ 中优化软件，它简化为数据在缓存中的访问和存储方式。

有关术语和详细信息，请参阅有关缓存的 wiki 条目，我将在此处缩小范围。

缓存按集和行进行组织。一次只使用一个集合，其中任何一个包含的行都可以使用。一行可以镜像的内存乘以行数，得出缓存大小。

对于一个特定的内存地址，我们可以用以下公式计算出哪个集合应该镜像它：

set = ( address / lineSize ) % numberOfsets

理想情况下，这种公式在集合之间给出了均匀的分布，因为每个内存地址都可能被读取（我说得很理想）。

很明显，可能会发生重叠。如果缓存未命中，则在缓存中读取内存并替换旧值。请记住，每组都有若干行，其中最近使用最少的行将被新读取的内存覆盖。

我将尝试在某种程度上遵循 Agner 的例子：

假设每个集合有 4 行，每行包含 64 个字节。我们首先尝试读取地址，该地址在集合中。然后我们还尝试读取地址 、 和 。所有这些都属于同一组。在阅读之前，集合中的所有行都会被占用。读取该内存会逐出集合中的现有行，即最初持有的行。问题在于，我们读取的地址（在本例中）是分开的。这是关键的一步（同样，对于这个例子）。0x2710280x2F000x37000x3F000x47000x47000x27100x800

临界步幅也可以计算：

criticalStride = numberOfSets * lineSize

间隔或多个相隔的变量争用相同的缓存行。criticalStride

这是理论部分。接下来是解释（也是Agner，我正在密切关注它以避免犯错误）：

假设一个 64x64 的矩阵（请记住，效果因缓存而异），缓存为 8kb，每组 4 行 * 行大小为 64 字节。每行可以容纳矩阵中的 8 个元素（64 位）。int

临界步幅为 2048 字节，对应于矩阵的 4 行（在内存中是连续的）。

假设我们正在处理第 28 行。我们尝试获取此行的元素，并将它们与第 28 列中的元素交换。该行的前 8 个元素组成一个缓存行，但它们将进入第 28 列中的 8 个不同的缓存行。请记住，临界步幅相距 4 行（一列中有 4 个连续元素）。

当列中达到元素 16 时（每组 4 行缓存行，相隔 4 行 = 麻烦），ex-0 元素将从缓存中逐出。当我们到达列的末尾时，所有先前的缓存行都将丢失，并且需要在访问下一个元素时重新加载（整行被覆盖）。

如果大小不是临界步幅的倍数，就会打乱这个完美的灾难场景，因为我们不再处理在垂直方向上相距临界步幅的元素，因此缓存重新加载的次数会大大减少。

另一个免责声明 - 我只是对解释有所了解，希望我搞定了，但我可能弄错了。无论如何，我正在等待 Mysticial 的回复（或确认）。:)

Luchian 解释了为什么会发生这种行为，但我认为展示这个问题的一个可能的解决方案，同时展示一些关于缓存遗忘算法的信息是个好主意。

你的算法基本上可以做到：

for (int i = 0; i < N; i++) 
   for (int j = 0; j < N; j++) 
        A[j][i] = A[i][j];

这对于现代 CPU 来说太可怕了。一种解决方案是了解有关缓存系统的详细信息并调整算法以避免这些问题。只要你知道这些细节，效果很好。不是特别便携。

我们能做得更好吗？是的，我们可以：解决这个问题的一般方法是缓存遗忘算法，顾名思义，它避免了依赖于特定的缓存大小 [1]

解决方案如下所示：

void recursiveTranspose(int i0, int i1, int j0, int j1) {
    int di = i1 - i0, dj = j1 - j0;
    const int LEAFSIZE = 32; // well ok caching still affects this one here
    if (di >= dj && di > LEAFSIZE) {
        int im = (i0 + i1) / 2;
        recursiveTranspose(i0, im, j0, j1);
        recursiveTranspose(im, i1, j0, j1);
    } else if (dj > LEAFSIZE) {
        int jm = (j0 + j1) / 2;
        recursiveTranspose(i0, i1, j0, jm);
        recursiveTranspose(i0, i1, jm, j1);
    } else {
    for (int i = i0; i < i1; i++ )
        for (int j = j0; j < j1; j++ )
            mat[j][i] = mat[i][j];
    }
}

稍微复杂一些，但一个简短的测试显示了我古老的 e8400 和 VS2010 x64 版本上非常有趣的东西，测试代码MATSIZE 8192

int main() {
    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);
    QueryPerformanceCounter(&start);
    recursiveTranspose(0, MATSIZE, 0, MATSIZE);
    QueryPerformanceCounter(&end);
    printf("recursive: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));

    QueryPerformanceCounter(&start);
    transpose();
    QueryPerformanceCounter(&end);
    printf("iterative: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));
    return 0;
}

results: 
recursive: 480.58ms
iterative: 3678.46ms

编辑：关于大小的影响：虽然在某种程度上仍然很明显，但它不那么明显，这是因为我们使用迭代解决方案作为叶节点，而不是递归到 1（递归算法的通常优化）。如果我们设置 LEAFSIZE = 1，缓存对我没有影响 [ - 在误差范围内，波动在 100 毫秒范围内;如果我们想要完全准确的值，这个“基准”我不会太舒服]）8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms

[1] 这些东西的来源：好吧，如果你不能从与 Leiserson 等人合作过的人那里得到讲座......我认为他们的论文是一个很好的起点。这些算法仍然很少被描述——CLR 有一个关于它们的脚注。尽管如此，这仍然是给人们带来惊喜的好方法。

---

编辑（注意：我不是发布此答案的人;我只是想添加这个）：
这是上述代码的完整 C++ 版本：

template<class InIt, class OutIt>
void transpose(InIt const input, OutIt const output,
    size_t const rows, size_t const columns,
    size_t const r1 = 0, size_t const c1 = 0,
    size_t r2 = ~(size_t) 0, size_t c2 = ~(size_t) 0,
    size_t const leaf = 0x20)
{
    if (!~c2) { c2 = columns - c1; }
    if (!~r2) { r2 = rows - r1; }
    size_t const di = r2 - r1, dj = c2 - c1;
    if (di >= dj && di > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, (r1 + r2) / 2, c2);
        transpose(input, output, rows, columns, (r1 + r2) / 2, c1, r2, c2);
    }
    else if (dj > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, r2, (c1 + c2) / 2);
        transpose(input, output, rows, columns, r1, (c1 + c2) / 2, r2, c2);
    }
    else
    {
        for (ptrdiff_t i1 = (ptrdiff_t) r1, i2 = (ptrdiff_t) (i1 * columns);
            i1 < (ptrdiff_t) r2; ++i1, i2 += (ptrdiff_t) columns)
        {
            for (ptrdiff_t j1 = (ptrdiff_t) c1, j2 = (ptrdiff_t) (j1 * rows);
                j1 < (ptrdiff_t) c2; ++j1, j2 += (ptrdiff_t) rows)
            {
                output[j2 + i1] = input[i2 + j1];
            }
        }
    }
}

为了说明Luchian Grigore的答案中的解释，以下是64x64和65x65矩阵两种情况下的矩阵缓存存在的样子（有关数字的详细信息，请参阅上面的链接）。

以下动画中的颜色表示以下含义：

• – 不在缓存中，
• – 在缓存中，
• – 缓存命中，
• – 只需从 RAM 中读取，
• – 缓存未命中。

64x64 案例：

请注意，几乎每次访问新行都会导致缓存未命中。现在它看起来是正常情况下的样子，一个 65x65 的矩阵：

在这里，您可以看到初始预热后的大多数访问都是缓存命中。这就是 CPU 缓存的一般工作方式。

---

_{可以在此处查看为上述动画生成帧的代码。}