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

Question

出于好奇,我决定将我自己的矩阵乘法函数与BLAS实现进行比较......我对结果的评价最少:

自定义实现,1000x1000矩阵乘法的10次试验:

Took: 15.76542 seconds.

BLAS实施,1000x1000矩阵乘法的10次试验:

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. 为什么我的实现速度要慢得多？

Answer 1

一个很好的出发点是Robert A. van de Geijn和Enrique S. Quintana-Ortí撰写的伟大的编程科学计算矩阵计算.它们提供免费下载版本.

BLAS分为三个级别:

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

• 2级函数是矩阵向量运算,例如一些矩阵向量乘积.这些功能可以用Level1功能实现.但是,如果您可以提供使用具有共享内存的某些多处理器体系结构的专用实现,则可以提高此功能的性能.

• 3级功能是矩阵矩阵产品之类的操作.您可以再次使用Level2函数实现它们.但Level3函数对O(N ^ 2)数据执行O(N ^ 3)运算.因此,如果您的平台具有缓存层次结构,那么如果您提供缓存优化/缓存友好的专用实现,则可以提高性能.书中很好地描述了这一点.Level3功能的主要推动力来自缓存优化.这种提升大大超过了并行性和其他硬件优化的第二次提升.

顺便说一下,大多数(甚至全部)高性能BLAS实现都没有在Fortran中实现.ATLAS是用C.实现的.GotoBLAS/OpenBLAS是用C语言实现的,它在Assembler中的性能关键部分.在Fortran中只实现了BLAS的参考实现.但是,所有这些BLAS实现都提供了Fortran接口,使其可以与LAPACK链接(LAPACK从BLAS获得所有性能).

优化的编译器在这方面起着次要作用(对于GotoBLAS/OpenBLAS,编译器根本不重要).

恕我直言,没有BLAS实现使用像Coppersmith-Winograd算法或Strassen算法这样的算法.我不确定原因,但这是我的猜测:

• 也许它不可能提供这些算法的缓存优化实现(即你会失去更多然后你会赢)
• 这些算法在数值上不稳定.由于BLAS是LAPACK的计算内核,因此这是不行的.

编辑/更新:

这个主题的新的和突破性的论文是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论文一起,很容易理解像英特尔MKL这样的库如何能够获得这样的性能.为什么使用行或列主存储并不重要!

最后的基准是在这里(我们称之为我们的项目ulmBLAS):

ulmBLAS,BLIS,MKL,openBLAS和Eigen的基准

另一个编辑/更新:

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

高性能LU分解

(这种LU分解例如由Matlab用于求解线性方程组.)

我希望有时间扩展教程来描述和演示如何在PLASMA中实现LU分解的高度可扩展的并行实现.

好的,这里你去:编码缓存优化的并行LU分解

PS:我也做了一些改善uBLAS性能的实验.实际上提升非常简单(是的,玩文字:))uBLAS的表现:

关于uBLAS的实验.

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

BLAZE的实验.

Answer 2

首先,BLAS只是一个约50个函数的接口.有许多竞争的接口实现.

首先,我会提到很大程度上不相关的事情:

• Fortran vs C,没什么区别
• 高级矩阵算法,如Strassen,实现不使用它们,因为它们在实践中没有帮助

大多数实现以或多或少的方式将每个操作分解为小维矩阵或向量运算.例如,大的1000x1000矩阵乘法可以分成50x50矩阵乘法的序列.

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

• SIMD风格的说明
• 指令级并行
• 缓存意识

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

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

(提示:这就是为什么如果您使用ATLAS,最好手动为您的特定机器构建和调整库,然后使用预先构建的库.)

Answer 3

首先,矩阵乘法的算法比你使用的算法更有效.

其次,您的CPU一次可以执行多于一条指令.

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

第三,你的代码远非最佳:

• 您正在使用原始指针,这意味着编译器必须假设它们可能是别名.您可以指定特定于编译器的关键字或标志,以告诉编译器它们不是别名.或者,您应该使用除原始指针之外的其他类型,以解决问题.
• 您通过对输入矩阵的每行/每列执行一次天真遍历来颠倒缓存.在继续下一个块之前,您可以使用阻塞在矩阵的较小块上执行尽可能多的工作,该块适合CPU缓存.
• 对于纯粹的数字任务,Fortran几乎是无与伦比的,C++需要大量的哄骗才能达到类似的速度.这是可以做到,并且有(通常使用表达式模板),这表明它的几个库,但它的不平凡,它不只是发生.

Answer 4

我不清楚BLAS的实现,但是对于矩阵乘法有更高效的算法比O(n3)复杂度更好.众所周知的是Strassen算法

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_ 的包装。