C中的优化矩阵乘法
我正在尝试比较矩阵乘法的不同方法。 第一个是常规方法：
do { for (j = 0; j < i; j++) { for (k = 0; k < i; k++) { suma = 0; for (l = 0; l < i; l++) suma += MatrixA[j][l]*MatrixB[l][k]; MatrixR[j][k] = suma; } } } c++; } while (c<iteraciones); 第二个包括首先转置矩阵B然后按行进行乘法运算：
int f, co; for (f = 0; f < i; f++) { for ( co = 0; co < i; co++) { MatrixB[f][co] = MatrixB[co][f]; } } c = 0; do { for (j = 0; j < i; j++) { for (k = 0; k < i; k++) { suma = 0; for (l = 0; l < i; l++) suma += MatrixA[j][l]*MatrixB[k][l]; MatrixR[j][k] = suma; } } } c++; } while (c<iteraciones);
第二种方法应该更快，因为我们正在访问连续的内存插槽，但我没有在性能上获得显着改进。 难道我做错了什么？
我可以发布完整的代码，但我认为不需要。
每个程序员应该了解的内存 （pdf链接）由Ulrich Drepper提供了很多关于内存效率的好主意，但特别是，他使用矩阵乘法作为了解内存和使用这些知识可以加快这一过程的一个例子。 请参阅他的论文中的附录A.1，并阅读6.2.1节。 文章中的表6.2显示，他可以将其运行时间从一个天真实现的1000×1000矩阵的时间变为10％。
当然，他的最终代码非常多毛，并且使用了许多系统特定的东西和编译时调整，但是，如果你真的需要速度，阅读那篇论文并阅读他的实现绝对是值得的。
实现这一目标可能非常重要。 对大型矩阵特别重要的一种优化是平铺乘法以将内容保留在高速缓存中。 我曾经测量过12倍的性能差异，但是我特意选择了一个消耗了我的缓存倍数的矩阵大小（大约’97，因此缓存很小）。
有很多关于这个主题的文献。 一个起点是：
http://en.wikipedia.org/wiki/Loop_tiling
有关更深入的研究，以下参考资料，尤其是Banerjee书籍，可能会有所帮助：
[Ban93] Banerjee，Utpal，Loop Transformations for Restructuring Compilers：the Foundations，Kluwer Academic Publishers，Norwell，MA，1993。
[Ban94] Banerjee，Utpal，Loop Parallelization，Kluwer Academic Publishers，Norwell，MA，1994。
[BGS93] Bacon，David F.，Susan L. Graham和Oliver Sharp，加利福尼亚大学伯克利分校计算机科学部高性能计算编译器转换，技术报告No UCB / CSD-93-781。
[LRW91] Lam，Monica S.，Edward E. Rothberg和Michael E Wolf。 阻塞算法的缓存性能和优化，在第四届编程语言的建筑支持国际会议上，于1991年4月在加利福尼亚州圣克拉拉举行，63-74。
[LW91] Lam，Monica S.和Michael E Wolf。 循环变换理论和最大化并行性的算法，在IEEE并行和分布式系统中，1991,2（4）：452-471。
[PW86] Padua，David A.和Michael J. Wolfe，超级计算机高级编译器优化，ACM通信，29（12）：1184-1201,1986。
[Wolfe89] Wolfe，Michael J.优化超级计算机超级编译器，麻省理工学院出版社，剑桥，MA，1989。
[Wolfe96] Wolfe，Michael J.，并行计算高性能编译器，Addison-Wesley，CA，1996。
注意：你的第二个实现中有一个BUG
for (f = 0; f < i; f++) { for (co = 0; co < i; co++) { MatrixB[f][co] = MatrixB[co][f]; } }
当你做f = 0时，c = 1
MatrixB[0][1] = MatrixB[1][0];
你覆盖了MatrixB[0][1]而失去了这个价值！ 当循环达到f = 1时，c = 0
MatrixB[1][0] = MatrixB[0][1];
复制的值与已存在的值相同。
如果矩阵不够大或者您不重复操作很多次，您将看不到明显的差异。
如果矩阵是，例如1,000×1,000，你将开始看到改进，但我会说，如果它低于100×100，你不应该担心它。
此外，任何“改进”可能是毫秒级，除非yoy要么使用非常大的矩阵要么重复操作数千次。
最后，如果您更换的计算机速度更快，差异会更小！
您获得多大的改进将取决于：
- 缓存的大小
- 缓存行的大小
- 缓存的关联度
对于小型矩阵和现代处理器，极有可能在第一次触摸后， MatrixA和MatrixB中的数据几乎完全保留在缓存中。
只需要你尝试一下（但这只会对大型矩阵产生影响）：从内部循环中的乘法逻辑中分离出你的加法逻辑，如下所示：
for (k = 0; k < i; k++) { int sums[i];//I know this size declaration is illegal in C. consider //this pseudo-code. for (l = 0; l < i; l++) sums[l] = MatrixA[j][l]*MatrixB[k][l]; int suma = 0; for(int s = 0; s < i; s++) suma += sums[s]; }
这是因为当你写suma时，你最终会停止你的管道。 当然，大部分都是在寄存器重命名等处理，但由于我对硬件的理解有限，如果我想从代码中挤出每一盎司的性能，我会这样做，因为现在你不需要停止管道等待写入suma。 由于乘法比加法更昂贵，你想让机器尽可能地与它并行化，因此保存加法的档位意味着你在加法循环中等待的时间少于在乘法循环中等待的时间。
这只是我的逻辑。 在该领域拥有更多知识的其他人可能不同意。
您是否可以发布一些数据，比较一系列矩阵大小的两种方法？ 可能是您的期望不切实际，而您的第二个版本更快但您尚未完成测量。
在测量执行时间时，不要忘记包括转置矩阵B的时间。
您可能想要尝试的其他方法是将代码的性能与BLAS库中的等效操作的性能进行比较。 这可能无法直接回答您的问题，但它可以让您更好地了解您对代码的期望。
两个N * N矩阵相乘的计算复杂度为O（N ^ 3）。 如果使用可能已被MATLAB采用的O（N ^ 2.73）算法，性能将得到显着提高。 如果您安装了MATLAB，请尝试将两个1024 * 1024矩阵相乘。 在我的计算机上，MATLAB以0.7s完成它，但像你这样的天真算法的C \ C ++实现需要20秒。 如果您真的关心性能，请参阅低复杂算法。 我听说存在O（N ^ 2.4）算法，但它需要一个非常大的矩阵，以便可以忽略其他操作。
你不应该写矩阵乘法。 你应该依赖外部库。 特别是你应该使用BLAS库中的GEMM例程。 GEMM通常提供以下优化
闭塞
高效矩阵乘法依赖于阻塞矩阵并执行几个较小的阻塞乘法。 理想情况下，选择每个块的大小以很好地适应缓存，从而大大提高性能 。
调音
理想的块大小取决于底层的内存层次结构（缓存有多大？）。 因此，应针对每台特定计算机调整和编译库。 除其他外，这是由ATLAS实施的BLAS 。
assembly级别优化
矩阵乘法很常见，开发人员会手动优化它。 特别是这在GotoBLAS完成。
异构（GPU）计算
Matrix Multiply非常FLOP /计算密集型，使其成为在GPU上运行的理想选择。 cuBLAS和MAGMA是很好的候选人。
简而言之，密集线性代数是一个研究得很好的主题。 人们致力于改进这些算法。 你应该使用他们的工作; 它会让他们开心。
没那么特别但更好：
c = 0; do { for (j = 0; j < i; j++) { for (k = 0; k < i; k++) { sum = 0; sum_ = 0; for (l = 0; l < i; l++) { MatrixB[j][k] = MatrixB[k][j]; sum += MatrixA[j][l]*MatrixB[k][l]; l++; MatrixB[j][k] = MatrixB[k][j]; sum_ += MatrixA[j][l]*MatrixB[k][l]; sum += sum_; } MatrixR[j][k] = sum; } } c++; } while (c
如果您正在处理小数字，那么您提到的改进可以忽略不计。 此外，性能将根据您运行的硬件而有所不同。 但如果你正在处理数百万的数字，那么它会起作用。 来到该程序，你可以粘贴你写的程序。
很老的问题，但是我的opengl项目的当前实现：
typedef float matN[N][N]; inline void matN_mul(matN dest, matN src1, matN src2) { unsigned int i; for(i = 0; i < N^2; i++) { unsigned int row = (int) i / 4, col = i % 4; dest[row][col] = src1[row][0] * src2[0][col] + src1[row][1] * src2[1][col] + .... src[row][N-1] * src3[N-1][col]; } }
其中N被矩阵的大小替换。 因此，如果您乘以4x4矩阵，那么您使用：
typedef float mat4[4][4]; inline void mat4_mul(mat4 dest, mat4 src1, mat4 src2) { unsigned int i; for(i = 0; i < 16; i++) { unsigned int row = (int) i / 4, col = i % 4; dest[row][col] = src1[row][0] * src2[0][col] + src1[row][1] * src2[1][col] + src1[row][2] * src2[2][col] + src1[row][3] * src2[3][col]; } }
这个函数主要是最小化循环，但模数可能很重要......在我的计算机上，这个函数比三重循环乘法函数快大约50％。
缺点：
-
需要大量代码（例如mat3 x mat3，mat5 x mat5 ......的不同function）
-
不规则乘法所需的调整（例如mat3 x mat4）.....
一般来说，转置B 应该比天真的实现快得多，但代价是浪费另一个NxN值的内存。 我只花了一周的时间来挖掘矩阵乘法优化，到目前为止，绝对不屈不挠的赢家是这样的：
for (int i = 0; i < N; i++) for (int k = 0; k < N; k++) for (int j = 0; j < N; j++) if (likely(k)) /* #define likely(x) __builtin_expect(!!(x), 1) */ C[i][j] += A[i][k] * B[k][j]; else C[i][j] = A[i][k] * B[k][j];
这甚至比早期评论中提到的Drepper方法更好，因为无论底层CPU的缓存属性如何，它都能以最佳方式工作。 诀窍在于重新排序循环，以便以行主顺序访问所有三个矩阵。