2012-07-10 349 views
168

在不同尺寸的矩阵上进行了一些实验之后,出现了一种模式。总体上,转置大小为2^n的矩阵比转置大小2^n+1中的一个要慢。对于n的小数值,差异不是很大。然而为什么转置512x512的矩阵要比转置513x513的矩阵慢得多?

大的差异发生在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毫秒

这究竟是为什么?

+7

您的代码看起来对我不友好。 – CodesInChaos 2012-07-10 13:02:26

+3

@CodeInChaos,它是。 – 2012-07-10 13:02:56

+7

这是几乎相同的问题,这个问题:http://stackoverflow.com/questions/7905760/matrix-multiplication-small-difference-in-matrix-size-large-difference-in-timi – Mysticial 2012-07-10 13:30:51

回答

157

这个解释来自于Optimizing software in C++中的Agner Fog,它减少了数据如何被访问和存储在缓存中。

有关条款和详细信息,请参阅wiki entry on caching,我将在此处缩小范围。

高速缓存组织在。一次只能使用一套,其中包含的任何一行都可以使用。一行可以镜像的行数乘以行数给我们缓存大小。

对于特定的存储器地址,就可以计算出该设定它应该与式进行镜像中:

set = (address/lineSize) % numberOfsets 

这类式是使整个组理想地均匀分布,因为每个存储器地址是作为很可能会被阅读(我说理想情况下)。

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

我将设法有所遵循从昂纳的例子:

假定每个组有4行,每行保持64个字节。我们首先尝试读取地址0x2710,该地址在集合28中。然后我们也尝试读取地址0x2F000x37000x3F00和​​。所有这些属于同一组。在阅读​​之前,该集合中的所有行将被占用。读取该内存会清除集合中现有的一行,该行最初持有0x2710。问题在于我们读取的地址是(此例)0x800。这是关键步幅(再次,对于这个例子)。

临界步幅也可以计算:

criticaStride = 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的回复(或确认)。 :)

+0

哦,下一次。只需通过[Lounge](http://chat.stackoverflow.com/rooms/10/loungec)直接ping我就可以了。我没有在SO上找到每个名称的实例。 :)我只通过定期的电子邮件通知看到了这一点。 – Mysticial 2012-07-10 13:39:40

+0

@Mysticial @Luchian Grigore我的一位朋友告诉我,他在'Ubuntu 11.04 i386'上运行的'Intel core i3' pc显示与* gcc 4.6 *几乎相同的性能。对于我的电脑'Intel Core 2 Duo',* gingw gcc4.4 *,在'windows 7(32)'上运行。当我编译这个带有* gcc 4.6 *的旧电脑'intel centrino'时,它显示出了很大的不同。 'Ubuntu 12.04 i386'。 – 2012-09-27 01:58:17

+0

另请注意,地址相差4096倍的内存访问会错误地依赖于Intel SnB系列CPU。 (即页面内的相同偏移量)。当某些操作是存储时,这可以降低吞吐量,负载和商店的混合。 – 2016-03-18 01:52:44

64

Luchian给出解释为什么这种行为发生,但我认为这会是一个不错的主意,以显示一种可能的解决了这个问题,并在同一时间出示了一下有关缓存忘却的算法。

你的算法基本上没有:

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

这太可怕了现代CPU。一个解决方案是知道你的缓存系统的细节,并调整算法以避免这些问题。只要你知道那些细节,工作很棒..不是特别便携。

我们可以做得更好吗?是的,我们可以:这个问题的一般方法是cache oblivious algorithms,作为它的名字说避免依赖于特定的缓存大小[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]; 
    } 
} 

稍微复杂一些,但一个简短的测试显示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 

的东西在我的古E8400与VS2010发布的x64挺有意思的,testcode编辑:关于大小的影响:它是那么明显,虽然仍明显在一定程度上,这是因为我们将迭代解决方案用作叶节点,而不是递归到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]; 
      } 
     } 
    } 
} 
+2

如果比较不同大小的矩阵之间的时间,而不是递归和迭代。在指定大小的矩阵上尝试递归解决方案。 – 2012-07-10 13:28:45

+0

@Luchian既然你已经解释了*为什么*他看到了这种行为,我认为在一般情况下为这个问题引入一个解决方案是相当有趣的。 – Voo 2012-07-10 13:32:52

+0

因为,我在质疑为什么一个更大的矩阵需要更短的时间来处理,而不是寻找更快的算法... – 2012-07-10 13:34:33

8

作为对Luchian Grigore's answer中解释的说明,以下是64x64和65x65矩阵这两种情况下的矩阵缓存存在情况(请参阅上面的链接,了解有关数字的详细信息)。下面

色彩的动画含义如下:

  • white - 不在缓存中,
  • light-green - 在高速缓存中,
  • bright green - 高速缓存命中,
  • orange - 从RAM刚读,
  • red - 缓存未命中。

64×64的情况下:

cache presence animation for 64x64 matrix

注意如何几乎每一个访问缓存未命中的一个新行的结果。现在怎么它看起来正常情况下,一个65x65矩阵:

cache presence animation for 65x65 matrix

在这里你可以看到,大部分的初始磨合后访问的高速缓存命中。这就是CPU缓存一般用于如何工作的方式。

+0

伟大的插图! – 2018-01-05 11:07:16