为什么转置512×512的matrix要比转置513×513的matrix慢得多？
在不同大小的matrix上进行了一些实验之后，出现了一个模式。 不变地， 转置大小为2^n的matrix比转置大小为2^n+1的matrix要慢 。 对于n小值，差异不是很大。
然而，在512的值上会出现很大的差异（至less对我而言）
免责声明：我知道这个函数实际上并没有将matrix转置，因为元素的双重交换，但是没有什么区别。
遵循代码：
#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毫秒
为什么发生这种情况？
这个解释来自于Agile Fog在C ++中优化软件 ，它减less了数据如何被访问和存储在caching中。
有关条款和详细信息，请参阅关于caching的wiki条目 ，我将在此处缩小范围。
caching按组和行组织。 一次只能使用一套，其中包含的任何一行都可以使用。 一行可以镜像的行数乘以行数给我们的caching大小。
对于一个特定的内存地址，我们可以用下面的公式来计算它应该映射到哪个集合：
set = ( address / lineSize ) % numberOfsets
这种公式给出了理想的均匀分布，因为每个内存地址都可能被读取（我理想地说）。
很明显，重叠可能发生。 在高速caching未命中的情况下，将在高速caching中读取内存，并replace旧值。 请记住，每一组都有许多行，最近最less使用的行将被新读取的内存覆盖。
我会尽量有点跟随Agner的例子：
假设每个集合有4行，每行都有64个字节。 我们首先尝试读取地址0x2710 ，它在第28集中。 然后我们也尝试读地址0x3F00和0x4700 。 所有这些都属于同一组。 在读取0x4700之前，集合中的所有行将被占用。 读取该内存会0x2710集合中现有的一行，该行最初保存的是0x2710 。 问题在于，我们读地址（这个例子）是0x800 。 这是关键的一步 （对于这个例子来说）。
关键的步幅也可以被计算出来：
criticaStride = numberOfSets * lineSize
variablescriticalStride或者多个分开来竞争相同的caching行。
这是理论部分。 接下来，解释（也是Agner，我正在密切关注以避免犯错误）：
假设一个64×64的matrix（记住，效果因caching而异），一个8kb的caching，每行4行* 64行的字节数。 每行可以包含matrix中的8个元素（64位int ）。
关键的步幅将是2048个字节，这对应于matrix的4行（在存储器中是连续的）。
假设我们正在处理第28行。我们试图获取这一行的元素，并将它们与第28列的元素进行交换。行的前8个元素组成了一个caching行，但是它们将进入8个不同的行请注意，关键步幅相差4行（一列中有4个连续的元素）。
当在列中到达元素16时（每组4个高速caching行和相距4行=故障），ex-0元素将从高速caching中被逐出。 当我们到达列的末尾时，所有先前的caching行将会丢失，并且在访问下一个元素（整个行被覆盖）时需要重新加载。
有一个不是关键跨度的倍数的大小会混淆这个完美的灾难场景 ，因为我们不再处理在垂直方向上非常关键的元素，所以高速caching重新加载的数量会大大减less。
另一个免责声明 – 我只是把我的头解释，希望我钉了它，但我可能是错的。 无论如何，我在等待Mysticial的回应（或确认）。 🙂
Luchian给出了为什么会发生这种行为的解释，但是我认为对这个问题展示一个可能的解决scheme并同时展示一些关于caching遗忘的algorithm是个不错的主意。
你的algorithm基本上做到：
for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) A[j][i] = A[i][j];
这对于现代CPU来说太恐怖了。 一个解决scheme是知道你的caching系统的细节，并调整algorithm，以避免这些问题。 只要你知道这些细节，工作很好..不是特别便携。
我们能做得更好吗？ 是的，我们可以：这个问题的一般方法是caching遗忘algorithm ，正如名称所述，避免依赖于特定的caching大小[1]
解决scheme将如下所示：
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]; } }
稍微复杂一点，但是一个简短的testing显示了我的古老的e8400与VS2010 x64版本相当有趣，testing代码为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
编辑：关于大小的影响：虽然仍然在一定程度上仍然显而易见，这是因为我们使用迭代解决scheme作为叶节点，而不是recursion到1（recursionalgorithm的通常优化）。 如果我们设置LEAFSIZE = 1，caching对我没有影响[ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms – 在误差范围内，波动在100ms区域内; 如果我们想要完全准确的价值，这个“基准”并不是我所愿意的）
[1]这个东西的来源：那么，如果你不能从与莱森森合作过的人那里得到一个演讲，那么我认为他们的论文是一个很好的起点。 这些algorithm还很less被描述 – CLR对它们有一个脚注。 尽pipe如此，这也是让人惊喜的好方法。
编辑 （注意：我不是发布这个答案的人，我只是想补充一点）：
这是上面代码的一个完整的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]; } } } }
