BLAS如何获得如此极致的性能？

Question

出于好奇，我决定对比一下自己的矩阵乘法函数和BLAS实现...结果让我大吃一惊：

Custom Implementation, 10 trials of 1000x1000 matrix multiplication:

Took: 15.76542 seconds.

BLAS Implementation, 10 trials of 1000x1000 matrix multiplication:

Took: 1.32432 seconds.

这是使用单精度浮点数。

我的实现：

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

我有两个问题：

1. 考虑到矩阵乘法，例如：nxm * mxn 需要 n*n*m 次乘法，在上述情况下需要 1000^3 或 1e9 次操作。在我的 2.6Ghz 处理器上，为什么 BLAS 可以在 1.32 秒内执行 10*1e9 次操作？即使乘法是单个操作且没有其他操作，也应该需要 ~4 秒。
2. 为什么我的实现速度如此之慢？

- DeusAduro

---

BLAS已经由该领域的专家进行了上下优化。我认为它正在利用您芯片上的SIMD浮点单元并采取许多技巧来改善缓存行为... - dmckee --- ex-moderator kitten

一个每秒运行2.63亿次的处理器如何在1.3秒内完成1千亿次操作？ - DeusAduro

多个执行单元、流水线技术以及单指令多数据（SIMD）技术可以同时对多组操作数进行相同的操作。一些编译器可以针对常见芯片上的SIMD单元进行优化，但您通常需要明确打开此功能，并且了解如何使用它（http://en.wikipedia.org/wiki/SIMD）。避免缓存未命中几乎肯定是最棘手的问题。 - dmckee --- ex-moderator kitten

假设是错误的。已知有更好的算法，请参见维基百科。 - MSalters

@DeusAduro：在我的回答中，针对如何编写一个可以与Eigen竞争的矩阵矩阵乘积？问题，我发布了一个小例子，展示了如何实现高效缓存的矩阵矩阵乘积。 - Michael Lehn

Answer 1

一个很好的起点是由Robert A. van de Geijn和Enrique S. Quintana-Ortí所写的伟大书籍The Science of Programming Matrix Computations。他们提供了免费下载版本。

BLAS被分为三个级别：

• Level 1 定义了一组仅对向量进行操作的线性代数函数。这些函数受益于矢量化（例如使用 SIMD，如 SSE）。

• Level 2 函数是矩阵-向量运算，例如某些矩阵-向量乘积。这些函数可以用 Level 1 函数来实现。但是，如果您能提供一个利用共享内存的多处理器架构的专用实现，则可以提高这些函数的性能。

• Level 3 函数是像矩阵-矩阵乘积这样的操作。同样，您可以用 Level 2 函数来实现它们。但是，Level 3 函数对 O(N^2) 数据执行 O(N^3) 操作。因此，如果您的平台具有缓存层次结构，则可以通过提供经过缓存优化/友好的专用实现来提高性能。这在书中有很好的描述。Level 3 函数的主要提升来自缓存优化。这种提升显著超过并行性和其他硬件优化的第二个提升。

顺便提一下，大多数（甚至全部）高性能BLAS实现都不是用Fortran编写的。ATLAS是用C实现的，GotoBLAS/OpenBLAS是用C实现的，其性能关键部分是用汇编语言实现的。只有BLAS的参考实现是用Fortran实现的。然而，所有这些BLAS实现都提供了Fortran接口，以便可以与LAPACK链接（LAPACK从BLAS获得所有性能）。

在这方面，优化编译器的作用很小（对于GotoBLAS/OpenBLAS，编译器根本无关紧要）。

在我看来，没有BLAS实现使用Coppersmith-Winograd算法或Strassen算法。可能的原因是：

• 也许不能提供这些算法的缓存优化实现（即你会失去更多而不是获得更多）
• 这些算法在数值上不稳定。由于BLAS是LAPACK的计算核心，所以这是不可接受的。
• 虽然这些算法在纸上具有不错的时间复杂度，但是大O表示法隐藏了一个很大的常数，因此只有对于非常大的矩阵才开始变得可行。

编辑/更新：

这个主题的新颖和开创性论文是BLIS论文。它们写得非常好。在我的“高性能计算软件基础”讲座中，我按照他们的论文实现了矩阵乘积。实际上，我实现了几个矩阵乘积的变体。最简单的变体完全使用普通的C语言编写，代码行数不到450行。所有其他变体仅仅优化循环。

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

矩阵乘法的整体性能仅取决于这些循环。大约99.9%的时间都在这里花费。在其他变体中，我使用了内部函数和汇编代码来提高性能。您可以在此处查看涵盖所有变体的教程：

ulmBLAS：GEMM（矩阵乘法）教程

结合BLIS论文，就可以很容易地理解像Intel MKL这样的库如何获得这样的性能。以及为什么使用行主还是列主存储不重要！

最终的基准测试结果在这里（我们称之为ulmBLAS）：

{{link2：ulmBLAS、BLIS、MKL、openBLAS和Eigen的基准测试}}

另一个编辑/更新：

我还撰写了一些关于如何使用BLAS解决数值线性代数问题（如解线性方程组）的教程：

高性能LU分解

（例如，Matlab使用此LU分解来解决线性方程组。）

我希望有时间扩展教程，描述并演示如何实现类似PLASMA的高度可扩展的并行LU分解。

好的，这里是：编写缓存优化的并行LU分解

P.S.：我还进行了一些改进uBLAS的实验。实际上，提高uBLAS的性能相当简单：

关于uBLAS的实验。

这是一个与BLAZE类似的项目：

在BLAZE上进行的实验。

- Michael Lehn

---

“Benchmarks for ulmBLAS, BLIS, MKL, openBLAS and Eigen” 的新链接：http://apfel.mathematik.uni-ulm.de/~lehn/ulmBLAS/#toc3 - Ahmed Fasih

原来IBM的ESSL使用了Strassen算法的变体 - https://www.ibm.com/support/knowledgecenter/zh/SSFHY8/essl_welcome.html - ben-albrecht

大多数链接都已失效 - Aurélien Pierre

TSoPMC的PDF可以在作者的页面上找到，网址为https://www.cs.utexas.edu/users/rvdg/tmp/TSoPMC.pdf。 - Alex Shpilkin

尽管Coppersmith-Winograd算法在理论上具有良好的时间复杂度，但大O符号隐藏了一个非常大的常数，因此它只对于非常大的矩阵才开始变得可行。 - Nihar Karve

Answer 2

首先，BLAS只是大约50个函数的接口。有许多竞争的接口实现。

首先我会提到一些与此无关的事情：

• Fortran和C之间无区别
• 高级矩阵算法，如Strassen，在实践中并没有帮助，因此实现不使用它们

大多数实现将每个操作分解为更小维度的矩阵或向量操作，这是比较明显的方式。例如，一个大的1000x1000矩阵乘法可能被分解成一系列50x50矩阵乘法。

这些固定大小的小维度操作（称为内核）使用针对其目标的多个CPU功能在特定CPU的汇编代码中进行硬编码：

• SIMD式指令
• 指令级并行性
• 缓存感知

此外，这些内核可以相对于彼此并行执行，使用多个线程（CPU内核），采用典型的map-reduce设计模式。

看一下ATLAS，它是最常用的开源BLAS实现。它有许多不同的竞争内核，并且在ATLAS库构建过程中，它在它们之间进行比赛（有些甚至是参数化的，因此相同内核可以具有不同的设置）。它尝试不同的配置，然后为特定目标系统选择最佳配置。

（提示：这就是为什么如果您使用ATLAS，则最好针对特定机器手动构建和调整库，而不是使用预构建的库。）

- Andrew Tomazos

---

ATLAS不再是最常用的开源BLAS实现。它已被OpenBLAS（GotoBLAS的分支）和BLIS（GotoBLAS的重构）超越。 - Robert van de Geijn

@ulaff.net：可能是这样。这篇文章是6年前写的。我认为目前最快的BLAS实现（当然是在Intel上）是Intel MKL，但它不是开源的。 - Andrew Tomazos

我同意你回答的精神。这里有一个学术链接，但它表明一些人已经使用Strassen类型/Winograd类型算法来实现真实世界的加速。https://www.ics.uci.edu/~paolo/FastMM/FMM-Reference/reference.html - creanion

Answer 3

首先，有比你正在使用的矩阵乘法更有效的算法。

其次，您的CPU可以同时处理多个指令。

您的CPU每个周期执行3-4条指令，并且如果使用了SIMD单元，则每个指令会处理4个浮点数或2个双精度浮点数。（当然，这个数字也不准确，因为CPU通常只能在一个周期内处理一个SIMD指令）

第三，您的代码远非最优：

• 您正在使用原始指针，这意味着编译器必须假设它们可能别名。您可以指定特定于编译器的关键字或标志来告诉编译器它们不别名。或者，您应该使用除原始指针之外的其他类型，这些类型会解决问题。
• 通过对输入矩阵的每行/列进行简单遍历，您正在使缓存失效。您可以使用分块将尽可能多的工作执行在矩阵的较小块上，该块适合CPU缓存，然后再转移到下一个块。
• 对于纯数值任务，Fortran几乎是无与伦比的，而要使C ++达到类似的速度需要大量努力。虽然这是可行的，并且有一些库可以证明它（通常使用表达式模板），但这并不是微不足道的，也不是“轻而易举”的。

- jalf

---

谢谢，根据Justicle的建议，我已经添加了限制正确代码的功能，但并没有看到很大的改进，我喜欢按块处理的想法。出于好奇，如果不知道CPU的缓存大小，怎样能写出优化的代码呢？ - DeusAduro

你不需要这样做。为了获得最佳代码，你需要知道CPU的缓存大小。 当然，这样做的缺点是你实际上在为一种CPU家族的最佳性能硬编码你的代码。 - jalf

至少在这里，内部循环避免了跨步加载。看起来这是为一个矩阵已经被转置而编写的。这就是为什么它比 BLAS "只"慢一个数量级！但是，是的，仍然会因缺乏高速缓存而崩溃。你确定 Fortran 会有很大帮助吗？我认为你所获得的只是 restrict(无别名) 是默认的，不像在 C / C++ 中。(不幸的是，ISO C++ 没有 restrict 关键字，在提供其作为扩展的编译器上，必须使用 __restrict__。) - Peter Cordes

Answer 4

我不了解具体的BLAS实现，但是有比O(n3)复杂度更好的矩阵乘法更高效的算法。其中一个众所周知的算法是Strassen Algorithm。

- softveda

---

Strassen算法在数值计算中没有被使用，原因有两个：1）它不稳定；2）虽然可以节省一些计算量，但代价是无法利用缓存层次结构，在实践中甚至会降低性能。 - Michael Lehn

针对于 Strassen 算法在 BLAS 库源码基础上的实际应用，最近有一篇论文发布：《Strassen Algorithm Reloaded》（链接：http://dl.acm.org/citation.cfm?id=3014983），在 SC16 会议上发表。该算法在问题规模为 1000x1000 时能够获得比 BLAS 更高的性能表现。 - Jianyu Huang

Answer 5

第二个问题的大部分论据，如汇编程序、分块等（但不包括低于N^3的算法，它们实在是过度开发了）都是有作用的。 但您的算法速度缓慢本质上是由矩阵大小和三个嵌套循环的不幸排列所致。 您的矩阵太大了，一次无法全部放入缓存内存中。 您可以重新安排循环，使尽可能多的操作在缓存中的一行上进行，从而显着减少缓存刷新（顺便说一下，将其分成小块具有类似的效果，最好的是相似地排列块上的循环）。 这里提供一个正方形矩阵的模型实现。 在我的计算机上，与标准实现（如您的实现）相比，它的时间消耗约为1：10。 换句话说：永远不要按我们在学校中学到的“行乘以列”的方案编写矩阵乘法。 在重新排列循环之后，通过展开循环、汇编代码等可以获得更多的改进。

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

还有一点需要注意的是：在我的电脑上，这个实现比替换所有内容为BLAS例程cblas_dgemm要更好（请在您的电脑上尝试！）。但是直接调用Fortran库中的dgemm_会快得多（1:4）。我认为这个例程实际上并不是Fortran而是汇编代码（我不知道库中有什么，我没有源代码）。对我来说完全不清楚为什么cblas_dgemm不像我所知道的那样快，因为它只是dgemm_的一个包装器。

- Wolfgang Jansen

---

比起我最初天真的基于数组的代码，速度快多了。太酷了！但是...即使我关闭了Julia的OpenMP多线程，仍然大约比Julia内置的BLAS矩阵乘法慢10倍。在我的M2 MacBook Air上测试了从1000x1000到5000x5000的矩阵大小。 - undefined

Answer 6

针对 MM 乘法中的原始代码，大多数操作的内存引用是性能不佳的主要原因。内存运行速度比缓存慢100-1000倍。

大多数加速来自于为此三重循环函数采用循环优化技术的应用。使用了两种主要的循环优化技术：展开和分块。关于展开，我们展开最外层的两个循环并将其分块以便在缓存中重复使用数据。外部循环展开通过减少整个操作期间对相同数据的不同时间的内存引用次数来优化数据访问时间。将循环索引阻塞在特定数字处有助于保留缓存中的数据。您可以选择优化 L2 缓存或 L3 缓存。

https://en.wikipedia.org/wiki/Loop_nest_optimization

- Pari Rajaram

Answer 7

这是一个实际的加速。关于使用SIMD汇编器优化C++代码的示例，请参见一些iPhone矩阵函数 - 这些函数比C版本快8倍以上，甚至没有进行“优化”的汇编 - 尚未进行流水线处理，并且存在不必要的堆栈操作。

此外，您的代码不符合"限制正确性" - 编译器如何知道修改C时，它不会修改A和B？

- Justicle

---

如果你像这样调用函数 mmult(A..., A..., A);，那么你肯定不会得到预期的结果。但我并不是在试图打败/重新实现BLAS，只是想看看它真正有多快，所以错误检查并没有考虑在内，只有基本功能。 - DeusAduro

抱歉，我想要清楚地表达的是，如果你在指针上加上“restrict”，你的代码会运行得更快。这是因为每次修改C时，编译器不必重新加载A和B，从而大大加快了内部循环的速度。如果你不相信，请查看反汇编结果。 - Justicle

@DeusAduro：这不是错误检查 - 可能编译器无法优化内部循环中对B[]数组的访问，因为它可能无法确定A和C指针从未与B数组重叠。如果存在别名，则在内部循环执行时B数组中的值可能会发生更改。将对B[]值的访问提升到内部循环之外并将其放入本地变量中，可能使编译器避免对B[]的持续访问。 - Michael Burr

谢谢你的澄清，我想我确实没有理解你的意思，我会看看能否添加并更新时间。 - DeusAduro

嗯，我首先尝试在VS 2008中使用“__restrict”关键字，应用于A、B和C。但结果没有改变。然而，将对B的访问从最内层循环移动到外部循环，时间提高了约10%。 - DeusAduro

抱歉，我对VC不确定，但是使用GCC时需要启用-fstrict-aliasing。这里也有更好的“restrict”解释：http://cellperformance.beyond3d.com/articles/2006/05/demystifying-the-restrict-keyword.html - Justicle

Answer 8

有很多原因。

首先，Fortran编译器高度优化，语言允许它们这样做。C和C++在数组处理方面非常松散（例如指针引用相同的内存区域），这意味着编译器无法预先知道该怎么做，被迫创建通用代码。在Fortran中，您的情况更加简化，编译器对发生的情况有更好的控制，使其能够进行更多的优化（例如使用寄存器）。

另一件事是Fortran按列存储数据，而C按行存储数据。我没有检查过您的代码，但要注意如何执行乘积。在C中，您必须按行扫描：这样可以沿着连续的内存扫描数组，减少缓存未命中。缓存未命中是效率低下的第一个来源。

第三，这取决于您正在使用的BLAS实现。某些实现可能是用汇编语言编写的，并针对您正在使用的特定处理器进行了优化。netlib版本是用Fortran 77编写的。

此外，您正在执行大量操作，其中大部分是重复和冗余的。所有这些乘法以获得索引都会损害性能。我不太清楚在BLAS中如何完成此操作，但有很多技巧可以防止昂贵的操作。

例如，您可以以这种方式重新编写您的代码。

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

试试吧，我相信你会节省一些东西。

关于你的第一个问题，原因是如果使用简单算法，矩阵乘法的复杂度为O(n^3)。但是有一些算法可以更好地扩展。

- Stefano Borini

---

抱歉，这个回答完全错误。BLAS实现不是用Fortran编写的。性能关键代码使用汇编语言编写，而现在最常见的是使用C语言进行编写。此外，BLAS将行/列顺序作为接口的一部分进行规定，实现可以处理任何组合。 - Andrew Tomazos

是的，这个回答完全是错误的。不幸的是，它充满了常见的无意义陈述，比如声称 BLAS 之所以更快是因为使用了 Fortran。拥有20个积极评价是一件坏事。现在，由于Stackoverflow的流行，这种无意义的说法甚至会进一步扩散！ - Michael Lehn

我认为你将未经优化的参考实现与生产实现混淆了。 参考实现只是用于指定库的接口和行为，并且由于历史原因是用Fortran编写的。 它不适用于生产使用。 在生产中，人们使用优化后的实现，这些实现表现出与参考实现相同的行为。 我已经研究了ATLAS（支持Octave-Linux“MATLAB”）的内部情况，我可以亲自确认它是用C / ASM编写的。商业实现几乎肯定也是如此。 - Andrew Tomazos

@KyleKanos：是的，这里是ATLAS的来源：http://sourceforge.net/projects/math-atlas/files/Stable/3.10.1/。据我所知，它是最常用的开源便携式BLAS实现。它是用C/ASM编写的。像英特尔这样的高性能CPU制造商也提供了专门针对其芯片进行优化的BLAS实现。我保证英特尔库的低级部分是用（duuh）x86汇编语言编写的，并且我相当确定中级部分将是用C或C++编写的。 - Andrew Tomazos

@KyleKanos：你有点困惑了。Netlib BLAS是参考实现，参考实现比优化实现要慢得多（请参见性能比较）。当有人说他们在使用一个集群上的netlib BLAS时，这并不意味着他们实际上在使用netlib参考实现。那样只是愚蠢的做法。这只是意味着他们正在使用一个与netlib blas具有相同接口的库。 - Andrew Tomazos