为什么矩阵乘法算法中的循环次序会影响性能?

我有两个函数来查找两个矩阵的乘积:

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倍)。任何想法为什么?

我相信你所看到的是计算机内存层次结构中引用局部性的影响。

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

大多数程序的编写方式是,如果内存中的单个字节被读入内存,程序稍后也会从该内存区域读取多个不同的值。 因此,这些高速缓存通常被设计为当您从内存中读取单个值时,该单个值周围的内存块(通常介于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] 。 这意味着在内循环的每次迭代中,您可能没有缓存未命中。

总的来说,这意味着代码的第二个版本不太可能在循环的每次迭代中都有缓存未命中,而第一个版本几乎肯定会。 因此,正如您所见,第二个循环可能比第一个循环更快。

有趣的是,许多编译器开始拥有原型支持,用于检测代码的第二个版本比第一个版本更快。 有些人会尝试自动重写代码以最大化并行性。 如果您有紫龙书的副本,第11章将讨论这些编译器的工作原理。

此外,您可以使用更复杂的循环进一步优化此循环的性能。 例如,一种称为阻塞的技术可以通过将数组拆分为可以在缓存中保存更长时间的子区域来显着提高性能,然后在这些块上使用多个操作来计算整体结果。

希望这可以帮助!

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

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

可能第二个必须在内存中跳过更多来访问数组元素。 它也可能是其他东西 – 您可以检查已编译的代码以查看实际发生的情况。

除了内存的位置,还有编译器优化。 向量和矩阵运算的关键是循环展开。

 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 [i] [k]在内存中连续访问
  • 内存访问和乘法可以在CPU中流水线化(几乎同时)。

如果n不是4或6或8的倍数怎么办? (或编译器决定将其展开的任何内容)编译器为您处理这个整理。 ;)

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

您可以做的另一件事就是multithreading增加性能。 这可以在4核CPU上将性能提高3倍。

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

第二个例子是c [i] [j]在每次迭代时都在变化,这使得它更难以优化。