一个很好的起点是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 实现。只有 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 论文，可以很容易地理解英特尔 MKL 等库如何获得这样的性能。以及为什么使用行或列主要存储并不重要！
最终的基准测试在这里（我们将项目称为 ulmBLAS）：
ulmBLAS、BLIS、MKL、openBLAS 和 Eigen 的基准
另一个编辑/更新：
我还写了一些关于 BLAS 如何用于数值线性代数问题的教程，例如求解线性方程组：
高性能 LU 分解
（例如，Matlab 使用此 LU 分解来求解线性方程组。）
我希望抽出时间来扩展本教程，以描述和演示如何像在PLASMA中那样实现高度可扩展的 LU 分解并行实现。
好的，你开始吧：编码缓存优化的并行 LU 分解
PS：我也做了一些提高uBLAS性能的实验。实际上，提高 uBLAS 的性能非常简单（是的，玩文字游戏 :)）：
uBLAS 上的实验。
这里有一个与BLAZE类似的项目：
对 BLAZE 的实验。