C中的优化矩阵乘法

Question

我正在尝试比较矩阵乘法的不同方法.第一个是常规方法:

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);

第二种方法应该更快,因为我们正在访问连续的内存插槽,但我没有在性能上获得显着改进.难道我做错了什么？

我可以发布完整的代码,但我认为不需要.

Answer 1

每个程序员应该了解的内存(pdf链接)由Ulrich Drepper提供了很多关于内存效率的好主意,但特别是,他使用矩阵乘法作为了解内存和使用这些知识可以加快这一过程的一个例子.请参阅他的论文中的附录A.1,并阅读6.2.1节.文章中的表6.2显示,对于1000x1000矩阵,他可以从一个天真的实现时间中获得10%的运行时间.

当然,他的最终代码非常多毛,并且使用了许多系统特定的东西和编译时调整,但是,如果你真的需要速度,阅读那篇论文并阅读他的实现绝对是值得的.

Answer 2

做到这一点非常重要.一种对大型矩阵特别重要的优化是平铺乘法以将内容保留在缓存中.我曾经测量过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.

Answer 3

注意:您的第二个实现中有一个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];

复制的值与已存在的值相同.

Answer 4

您不应该编写矩阵乘法。您应该依赖外部库。特别是您应该使用库GEMM中的例程BLAS。GEMM通常提供以下优化

阻塞

高效矩阵乘法依赖于分块矩阵并执行几个较小的分块乘法。理想情况下，选择每个块的大小以很好地适应缓存，从而大大提高性能。

调音

理想的块大小取决于底层内存层次结构（缓存有多大？）。因此，应该针对每台特定机器调整和编译库。除其他外，这是通过ATLAS实施BLAS.

装配级优化

矩阵乘法如此常见，以至于开发人员会手动优化它。特别是这是在GotoBLAS.

异构(GPU)计算

矩阵乘法是 FLOP/计算密集型的，使其成为在 GPU 上运行的理想选择。 cuBLAS并且MAGMA是这方面的良好候选人。

简而言之，稠密线性代数是一个经过深入研究的课题。人们毕生致力于这些算法的改进。你应该使用他们的作品；这会让他们高兴。

Answer 5

如果矩阵不够大或者您没有多次重复操作，您将不会看到明显的差异。

如果矩阵是 1,000x1,000，您将开始看到改进，但我想说，如果它低于 100x100，您不必担心。

此外，任何“改进”都可能是毫秒的量级，除非您正在处理非常大的矩阵或重复操作数千次。

最后，如果您将正在使用的计算机更换为更快的计算机，则差异会更小！