为什么转置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集中。 然后我们也尝试读地址0x3F000x4700 。 所有这些都属于同一组。 在读取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]; } } } }