BLAS是如何获得如此极端的性能的?

一个好的起点是这本伟大的书 The Science of Programming Matrix Computations 罗伯特·范德盖恩和恩里克·昆塔纳·奥特合著。他们提供免费下载版本。

BLAS分为三个级别:

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

• 3级功能是与矩阵乘积类似的操作。同样,您可以根据Level2函数来实现它们。但是Level3函数对O(N^2)数据执行O(N^3)操作。因此,如果您的平台具有缓存层次结构,那么如果您提供一个专门的实现 缓存优化/缓存友好

顺便说一句,大多数(甚至全部)高性能BLAS实现都不是用Fortran实现的。ATLAS是在C中实现的。GotoBLAS/OpenBLAS是在C中实现的,其性能关键部件是在汇编程序中实现的。只有BLAS的参考实现是用Fortran实现的。然而,所有这些BLAS实现都提供了一个Fortran接口,这样它就可以与LAPACK链接(LAPACK从BLAS获得所有性能)。

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

IMHO没有BLAS实现使用诸如铜匠维诺格拉德算法或斯特拉森算法之类的算法。可能的原因是:

• 也许不可能为这些算法提供缓存优化的实现(即,你会失去更多,然后你会赢)
• 这些算法在数值上不稳定。由于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因式分解的高度可扩展并行实现,如 PLASMA .

好的,给你: Coding a Cache Optimized Parallel LU Factorization

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

Experiments on uBLAS

这里有一个类似的项目 BLAZE

Experiments on BLAZE

首先,BLAS只是一个大约50个函数的接口。该接口有许多相互竞争的实现。

• 先进的矩阵算法,如Strassen,实现不使用它们,因为它们在实践中没有帮助

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

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

• SIMD样式说明
• 指令级并行
• 缓存感知

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

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

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

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

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

• 您使用的是原始指针,这意味着编译器必须假设它们可能是别名。您可以指定编译器特定的关键字或标志,以告知编译器它们不使用别名。或者,您应该使用原始指针以外的其他类型来解决问题。
• 对于纯粹的数值任务,FORTRAN几乎是不可替代的,C++需要大量的哄骗来达到类似的速度。这是可以做到的,并且有一些库演示了它(通常使用表达式模板),但它不是琐碎的,也不是 只是 发生

我不知道BLAS的具体实现,但是有比O(n3)复杂度更好的更有效的矩阵乘法算法。众所周知的是 Strassen Algorithm

第二个问题的大多数参数——汇编程序、分解成块等(但不是少于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úU dgemm没有那么快,因为据我所知,它只是dgemmú的包装。

iPhone matrix functions -它们比C版本快8倍多,甚至还没有“优化”组装-还没有管道衬里,也没有不必要的堆栈操作。

此外,您的代码不是“ restrict correct “-编译器如何知道当它修改C时,它不是在修改A和B?

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

大部分的加速来自于在MM乘法中对这个三环函数采用循环优化技术。采用了两种主回路优化技术;展开和阻塞。关于展开,我们展开最外层的两个循环,并将其阻塞,以便在缓存中重用数据。外部循环展开通过减少在整个操作过程中不同时间对相同数据的内存引用数量,帮助优化临时数据访问。在特定编号处阻塞循环索引,有助于将数据保留在缓存中。您可以选择优化二级缓存或三级缓存。

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

另一件事是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;
           }
     }
  }
} 

试试看,我相信你会省下一些钱。

关于你的#1问题,原因是如果你使用一个简单的算法,矩阵乘法的比例是O(n^3)。有些算法 scale much better .