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

Question

在对不同大小的方阵进行了一些实验之后，出现了一种模式。转置一个 size 的矩阵2^n2^n+1总是比转置一个 size 的矩阵要慢。对于较小的值n，差异不大。

然而，在 512 的值上会出现很大的差异。（至少对我来说）

免责声明：我知道由于元素的双重交换，该函数实际上并没有转置矩阵，但这没有区别。

遵循代码：

#define SAMPLES 1000
#define MATSIZE 512

#include <time.h>
#include <iostream>
int mat[MATSIZE][MATSIZE];

void transpose()
{
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
   {
       int aux = mat[i][j];
       mat[i][j] = mat[j][i];
       mat[j][i] = aux;
   }
}

int main()
{
   //initialize matrix
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
       mat[i][j] = i+j;

   int t = clock();
   for ( int i = 0 ; i < SAMPLES ; i++ )
       transpose();
   int elapsed = clock() - t;

   std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed / SAMPLES;
}

改变MATSIZE让我们改变大小（呃！）。我在ideone上发布了两个版本：

大小 512 - 平均2.46 毫秒- http://ideone.com/1PV7m
大小 513 - 平均0.75 毫秒- http://ideone.com/NShpo

在我的环境中（MSVS 2010，完全优化），区别是相似的：

大小 512 - 平均2.19 毫秒
大小 513 - 平均0.57 毫秒

为什么会这样？

score 208 · Accepted Answer

解释来自 Agner Fog 在Optimizing software in C++ 中，它简化了数据在缓存中的访问和存储方式。

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

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

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

set = ( address / lineSize ) % numberOfsets

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

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

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

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

还可以计算临界步幅：

criticalStride = numberOfSets * lineSize

criticalStride间隔或多个分开的变量竞争相同的高速缓存行。

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

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

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

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

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

大小不是临界步长的倍数会使这个完美的灾难场景变得混乱，因为我们不再处理在垂直方向上临界步距分开的元素，因此缓存重新加载的次数大大减少。

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

score 85 · Accepted Answer

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

以下动画中的颜色含义如下：

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

64x64 案例：

请注意，几乎每次访问新行都会导致缓存未命中。现在它是如何寻找正常情况的，一个 65x65 矩阵：

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

_{为上述动画生成帧的代码可以在这里看到。}

score 80 · Accepted Answer

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，缓存对我没有影响 [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms- 在误差范围内，波动在 100ms 范围内；如果我们想要完全准确的值，我不会对这个“基准”感到满意]）

[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];
            }
        }
    }
}

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

3 回答 3

Related

Reference