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

在对不同大小的方阵进行了一些实验后，形成了一个模式。总是,将大小为__ABC0的矩阵转置比将大小为2^n+1的矩阵转置要慢。对于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上发布了两个版本:

在我的环境中(MSVS 2010，完全优化)，差异是类似的:

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

为什么会这样?

32692 次浏览
小开
最佳答案

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

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

缓存被组织在集和行中。一次只使用一个集合，其中包含的任何行都可以使用。一行能镜像的内存乘以行数就得到了缓存大小。

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

set = ( address / lineSize ) % numberOfsets

理想情况下，这种公式在集合上给出了统一的分布，因为每个内存地址都可能被读取(我说的是在理想的情况下)。

很明显，重叠是可能发生的。在缓存丢失的情况下，在缓存中读取内存并替换旧的值。请记住，每个集合都有一些行，其中最近最少使用的行将被新读取的内存覆盖。

我将试着遵循Agner的例子:

假设每个集合有4行，每行包含64字节。我们首先尝试读取地址0x2710，它位于集合28中。然后我们还尝试读取地址0x2F00， 0x3700， 0x3F00和0x4700。所有这些都属于同一个集合。在读取0x4700之前，集合中的所有行都将被占用。读取该内存将清除集合中现有的一行，即最初保存0x2710的行。问题在于我们读取的地址(在本例中)是0x800分开的。这是关键的一步(同样，对于这个例子)。

关键步幅也可以计算:

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版本，testcode for 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和co一起工作的人那里得到一个讲座。我认为他们的论文是一个很好的起点。这些算法仍然很少被描述——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矩阵两种情况下的矩阵缓存存在情况(有关数字的详细信息，请参阅上面的链接)。

下面动画中的颜色意味着:

  • white -不在缓存中，
  • 浅绿色 -在缓存中，
  • 鲜绿 -缓存命中，
  • orange -从RAM中读取，
  • red缓存丢失。

64x64机箱:

64x64矩阵的缓存存在动画

注意几乎每一个对新行的访问是如何导致缓存失败的。现在看看正常情况下65x65矩阵的情况:

缓存存在动画65x65矩阵

在这里，您可以看到在初始热身之后的大多数访问都是缓存命中。这就是CPU缓存的工作原理。

为以上动画生成帧的代码可以看到here。


提交答案