2011-09-13 30 views
15

我给定两个函数用于查找两个矩阵的乘积:为什么矩阵乘法算法中的循环顺序会影响性能?

void MultiplyMatrices_1(int **a, int **b, int **c, int n){ 
     for (int i = 0; i < n; i++) 
      for (int j = 0; j < n; j++) 
       for (int k = 0; k < n; k++) 
        c[i][j] = c[i][j] + a[i][k]*b[k][j]; 
    } 

void MultiplyMatrices_2(int **a, int **b, int **c, int n){ 
     for (int i = 0; i < n; i++) 
      for (int k = 0; k < n; k++) 
       for (int j = 0; j < n; j++) 
        c[i][j] = c[i][j] + a[i][k]*b[k][j]; 
} 

我跑和使用gprof,每个具有除了该功能相同的代码异形两个可执行。对于2048 x 2048大小的矩阵,其中第二个显着(约5倍)。关于为什么?

回答

30

我相信你在看的是locality of reference在计算机内存层次结构中的影响。

典型地,计算机存储器被分成具有不同的性能特性的不同类型(这通常被称为memory hierarchy)。最快的存储器位于处理器的寄存器中,可以(通常)在一个时钟周期内访问和读取。但是,这些寄存器通常只有少数(通常不超过1KB)。另一方面,计算机的主存很大(如8GB),但访问速度要慢得多。为了提高性能,计算机通常在物理上构造成在处理器和主存储器之间具有several levels of caches。这些缓存比寄存器慢,但比主内存快得多,所以如果你在缓存中查看内存访问,它往往比如果你必须去主内存快得多(通常在5-25倍之间更快)。当访问内存时,处理器首先检查内存缓存中的该值,然后返回主内存读取该值。如果始终访问缓存中的值,则最终性能会比跳过内存时的性能好得多内存,随机访问值。

大多数程序都是以一种方式写入的,如果内存中的单个字节被读入内存,程序稍后会从该内存区域读取多个不同的值。因此,这些缓存通常被设计成当你从内存中读取单个值时,围绕该单值的一块内存(通常介于1KB和1MB之间)的值也被拉入缓存。这样,如果你的程序读取附近的值,它们已经在缓存中,你不必去主内存。

现在,最后一个细节 - 在C/C++中,数组以行优先顺序存储,这意味着矩阵的单个行中的所有值都相邻存储。因此,在内存中,阵列看起来像第一行,然后是第二行,然后是第三行等。

鉴于此,我们来看看您的代码。第一个版本如下所示:

for (int i = 0; i < n; i++) 
     for (int j = 0; j < n; j++) 
      for (int k = 0; k < n; k++) 
       c[i][j] = c[i][j] + a[i][k]*b[k][j]; 

现在,我们来看看最内层的代码行。在每次迭代中,k的值都在变化。这意味着在运行最内层循环时,循环的每次迭代很可能在加载b[k][j]的值时发生缓存缺失。原因在于,因为矩阵以行优先顺序存储,所以每次增加k时,都会跳过矩阵的整行并进一步跳到内存中,可能远远超过已缓存的值。但是,在查找c[i][j]时(因为ij是相同的),您不会错过,也不会错过a[i][k],因为这些值是按行优先顺序排列的,并且如果a[i][k]的值是从先前缓存的迭代,在此迭代中读取的a[i][k]的值是来自相邻的存储器位置。因此,在最内层循环的每次迭代中,您可能会有一次高速缓存未命中。

但考虑这第二个版本:

for (int i = 0; i < n; i++) 
     for (int k = 0; k < n; k++) 
      for (int j = 0; j < n; j++) 
       c[i][j] = c[i][j] + a[i][k]*b[k][j]; 

现在,因为你在每次迭代增加j,让我们想想有多少高速缓存未命中你可能对最里面的语句。因为这些值是按行优先顺序排列的,所以c[i][j]的值很可能是缓存中的,因为前一次迭代中的c[i][j]的值可能会被缓存并准备好读取。同样,b[k][j]可能被缓存,并且由于ik不会更改,所以缓存a[i][k]的几率也是可能的。这意味着在内部循环的每次迭代中,您可能没有缓存未命中。总的来说,这意味着代码的第二个版本不太可能在循环的每次迭代中都存在缓存未命中,而第一个版本几乎肯定会发生。因此,如你所见,第二个循环可能比第一个循环更快。有趣的是,许多编译器开始有原型支持来检测代码的第二版本比第一版本更快。有些人会尝试自动重写代码以最大化并行性。如果您有Purple Dragon Book的副本,第11章将讨论这些编译器的工作方式。

此外,您可以使用更复杂的循环来进一步优化此循环的性能。例如,一种名为blocking的技术可用于通过将阵列拆分为可保存在缓存中的更长的子区域,然后对这些区块使用多个操作来计算总体结果,从而显着提高性能。

希望这会有所帮助!

+1

+1确实不错!另外,关于缓存调试的建议@Kerrek SB增加了更多的技术细节。 – rbaleksandar

0

可能是第二个不得不在内存中跳过访问数组元素。它也可能是其他的东西 - 你可以检查编译后的代码,看看实际发生了什么。

5

这可能是记忆的地方。当您重新排序循环时,最内层循环中所需的内存越来越近并且可以被缓存,而在低效版本中,您需要从整个数据集中访问内存。

测试此假设的方法是在两段代码上运行缓存调试器(如cachegrind),并查看它们产生多少缓存未命中。

+0

+1 cachegrind –

0

除了内存的地方之外,还有编译器优化。矢量和矩阵运算的关键之一是循环展开。

for (int k = 0; k < n; k++) 
    c[i][j] = c[i][j] + a[i][k]*b[k][j]; 

您可以在此内环ij不改看。这意味着它可以被重写为

for (int k = 0; k < n; k+=4) { 
    int * aik = &a[i][k]; 
    c[i][j] += 
     + aik[0]*b[k][j] 
     + aik[1]*b[k+1][j] 
     + aik[2]*b[k+2][j] 
     + aik[3]*b[k+3][j]; 
} 

你可以看到会有

  • 四次循环更少并且访问到c [i] [j]
  • A [1] [k]的正在内存中连续访问
  • 内存访问和倍增可以在CPU中流水线(几乎同时)。

如果n不是4或6或8的倍数? (或者编译器决定将其展开到哪里)编译器会为你处理这个问题。;)

为了加快此解决方案的速度,您可以尝试先转置b矩阵。这是一些额外的工作和编码,但这意味着对b转置的访问在内存中也是连续的。 (当你用[j]交换[k]时)

你可以做的另一件事是提高性能,以便多线程乘法。这可以在4核CPU上将性能提高3倍。

最后,你可能会考虑使用floatdouble你可能会认为int会更快,但是这是情况并非总是如此的浮点运算可以更加高度优化的(无论是在硬件和编译器)

第二例如,c [i] [j]在每次迭代中都会发生变化,这使得优化变得更加困难。