c++ - 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];
}

我有两个问题:

鉴于矩阵-矩阵乘法表示:nxm * mxn 需要 n*n*m 次乘法，因此在 1000^3 或 1e9 运算以上的情况下。我的 2.6Ghz 处理器如何让 BLAS 在 1.32 秒内完成 10*1e9 次操作？即使乘法是一个单一的操作并且没有做任何其他事情，它也应该需要大约 4 秒。

为什么我的实现速度这么慢？

Answer 1

好的起点是好书 The Science of Programming Matrix Computations作者:Robert A. van de Geijn 和 Enrique S. Quintana-Ortí。他们提供免费下载版本。
BLAS分为三个层次:

级别 1 定义了一组仅对 vector 进行运算的线性代数函数。这些函数受益于矢量化(例如使用 SSE)。

2 级函数是矩阵 vector 运算，例如一些矩阵 vector 乘积。这些功能可以按照Level1 功能来实现。但是，如果您可以提供一个专用实现来利用某些具有共享内存的多处理器体系结构，则可以提高此函数的性能。

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

顺便说一句，大多数(甚至全部)高性能 BLAS 实现都没有在 Fortran 中实现。 ATLAS 是用 C 实现的。 GotoBLAS/OpenBLAS 是用 C 实现的，它的性能关键部分在 Assembler 中实现。 Fortran中只实现了BLAS的引用实现。然而，所有这些 BLAS 实现都提供了一个 Fortran 接口(interface)，以便它可以链接到 LAPACK(LAPACK 从 BLAS 中获得所有性能)。
优化的编译器在这方面的作用很小(对于 GotoBLAS/OpenBLAS，编译器根本不重要)。
恕我直言，没有 BLAS 实现使用像 Coppersmith-Winograd 算法或 Strassen 算法这样的算法。可能的原因是:

也许不可能提供这些算法的缓存优化实现(即你会失去更多然后你会赢)

这些算法在数值上不稳定。由于 BLAS 是 LAPACK 的计算内核，因此这是行不通的。

尽管这些算法在纸面上具有很好的时间复杂度，但大 O 表示法隐藏了一个大常数，因此它只开始适用于非常大的矩阵。

编辑/更新:
该主题的新论文和开创性论文是 BLIS papers .他们写得特别好。在我的讲座“高性能计算的软件基础”中，我按照他们的论文实现了矩阵-矩阵产品。实际上我实现了矩阵矩阵产品的几个变体。最简单的变体完全用纯 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: Tutorial on GEMM (Matrix-Matrix Product)
与 BLIS 论文一起，很容易理解像英特尔 MKL 这样的库如何获得这样的性能。以及为什么使用行或列主要存储并不重要!
最终的基准测试在这里(我们称我们的项目为 ulmBLAS):
Benchmarks for ulmBLAS, BLIS, MKL, openBLAS and Eigen
另一个编辑/更新:
我还写了一些关于如何将 BLAS 用于求解线性方程组等数值线性代数问题的教程:
High Performance LU Factorization
(例如，此 LU 分解被 Matlab 用于求解线性方程组。)
我希望能抽出时间来扩展本教程，以描述和演示如何实现 LU 分解的高度可扩展并行实现，如 PLASMA .
好的，给你: Coding a Cache Optimized Parallel LU Factorization
P.S.:我也做了一些关于提高 uBLAS 性能的实验。提高 uBLAS 的性能实际上非常简单(是的，玩文字 :) ):
Experiments on uBLAS .
这里有一个类似的项目 BLAZE :
Experiments on BLAZE .