BLAS 是如何获得如此极端的表现的？

出于好奇，我决定将我自己的矩阵乘法功能与 BLAS 的实现进行比较... ... 我对结果表示最不惊讶的是:

定制实施，10个试验 1000x1000矩阵乘法:

Took: 15.76542 seconds.

推行《基本法》、10项试验计划 1000x1000矩阵乘法:

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.6 Ghz 处理器上，BLAS 怎么可能在1.32秒内完成10 * 1e9操作？即使乘法是一个单独的操作，并且没有其他操作，也需要大约4秒钟。
  2. 为什么我的实现要慢得多？
51536 次浏览
小开

因为很多原因。

首先，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)。有一些算法。

小开

这是一个现实的加速。举一个 SIMD 汇编程序在 C + + 代码上可以做什么的例子，看一些例子 IPhone 矩阵函数-这些比 C 版本快8倍多，甚至没有“优化”汇编-还没有管道衬里和不必要的堆栈操作。

而且你的代码不是“ 限制正确”——编译器怎么知道它修改 C 的时候不是在修改 A 和 B 呢？

小开

我不是很清楚 BLAS 的具体实现，但是有比 O (n3)复杂度更高效的矩阵乘法算法。最著名的是 施特拉森演算法

小开

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

其次，您的 CPU 一次可以执行多条指令。

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

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

  • 您使用的是原始指针，这意味着编译器必须假定它们可能是别名。可以指定特定于编译器的关键字或标志来告诉编译器它们不是别名。或者，您应该使用原始指针以外的其他类型，这些类型可以解决问题。
  • 通过对输入矩阵的每一行/每一列执行初始遍历，可以颠覆缓存。在移动到下一个块之前，您可以使用阻塞在适合 CPU 缓存的较小的矩阵块上执行尽可能多的工作。
  • 对于纯粹的数字任务，Fortran 几乎是无敌的，而 C + + 也需要很多努力才能达到类似的速度。这是可以做到的，并且有一些库对此进行了演示(通常使用表达式模板) ，但是这并不简单，而且不会发生 只是。
小开
最佳答案

一个很好的出发点是罗伯特 · A · 范德盖因和恩里克 · S · 昆塔纳-奥尔蒂的伟大著作《 编程矩阵计算科学》。他们提供免费下载版本。

营商联络小组分为三个层次:

  • 第1级定义了一组仅对向量进行操作的线性代数函数。这些函数受益于向量化(例如使用 SSE)。

  • 二级函数是矩阵向量运算，例如一些矩阵向量乘积。这些函数可以用 Level 1函数来实现。但是，如果可以提供一个专用的实现，使用一些具有共享内存的多处理器体系结构，则可以提高这些函数的性能。

  • 第3级函数是类似矩阵矩阵乘积的运算。同样，您可以根据 Level 2函数来实现它们。但 Level 3函数对 O (N ^ 2)数据执行 O (N ^ 3)操作。因此，如果您的平台有一个缓存层次结构，那么如果您提供一个专用的实现，即 缓存优化/缓存友好，那么您可以提高性能。这在书中有很好的描述。Level 3函数的主要增强来自缓存优化。这个增强大大超过了并行性和其他硬件优化带来的第二个增强。

顺便提一下，大部分(甚至全部)高性能的《基本法》实施计划并没有在 Fortran 实施。ATLAS 在 C 语言中实现 GotoBLAS/OpenBLAS 在 C 语言中实现，其性能关键部分在汇编语言中实现。Fortran 只实施《基本法》的参考实施办法。然而，所有这些 BLAS 实现都提供了一个 Fortran 接口，这样就可以将其链接到 LAPACK (LAPACK 通过 BLAS 获得所有性能)。

优化的编译器在这方面只起到很小的作用(对于 GotoBLAS/OpenBLAS，编译器根本不重要)。

恕我直言，没有 BLAS 的实现使用了像 Coppersmith-Winograd 算法或者施特拉森演算法这样的算法。可能的原因是:

  • 也许不可能提供这些算法的缓存优化实现(也就是说，你会失去比你会赢得更多)
  • 这些算法在数值上是不稳定的，因为 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) :

UlmBLAS、 BLIS、 MKL、 openBLAS 和 Eigen 的基准

另一个编辑/更新:

我还写了一些关于 BLAS 如何用于解决数值线性代数线性方程组等问题的教程:

高性能 LU 因子分解

(这个逻辑单位分解就是 Matlab 用来求解线性方程组的例子。)

我希望找到时间 来扩展本教程来描述和演示如何像 血浆那样实现一个高度可伸缩的 LU 分解并行实现。

好的，给你: 一种 Cache 优化并行 LU 分解编码算法

附注: 我也做了一些提高 uBLAS 性能的实验。实际上，提高 uBLAS 的性能非常简单(是的，玩文字游戏:) :

UBLAS 实验。

这里有一个与 BLAZE类似的项目:

BLAZE 实验。

小开

因此，首先 BLAS 只是一个约50个功能的接口。接口有许多相互竞争的实现。

首先，我要提到一些基本上不相关的事情:

  • Fortran 和 C 没有区别
  • 先进的矩阵算法，如施特拉森，实现不使用它们，因为它们在实践中没有帮助

大多数实现都会以或多或少明显的方式将每个操作分解为小维矩阵或向量操作。例如，一个1000x1000的大矩阵乘法可能会被分解成一个50x50的矩阵乘法序列。

这些固定大小的小维操作(称为内核)使用其目标的几个 CPU 特性在特定于 CPU 的汇编代码中进行硬编码:

  • SIMD 风格的指令
  • 指令层级平行
  • 缓存意识

此外，在典型的 map-reduce 设计模式中，这些内核可以使用多个线程(CPU 内核)相互并行地执行。

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

(提示: 这就是为什么如果您正在使用 ATLAS，您最好手工为您的特定机器构建和调优库，而不是使用预构建的库。)

小开

第二个问题的大多数参数——汇编程序、分块等(但不是小于 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 的一个包装器。

小开

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

大部分的加速来自采用环路优化技术的这三个环路函数在 MM 乘。主要采用两种循环优化技术: 展开和阻塞。关于展开，我们展开最外面的两个循环并阻塞它，以便在缓存中进行数据重用。外部循环展开通过在整个操作过程中的不同时间减少对相同数据的内存引用数量，从而帮助在时间上优化数据访问。将循环索引阻塞在特定的数字上，有助于将数据保留在缓存中。您可以选择优化 L2缓存或 L3缓存。

Https://en.wikipedia.org/wiki/loop_nest_optimization


提交答案