BLAS如何获得如此极致的performance？
出于好奇，我决定将自己的matrix乘法函数与BLAS实现进行比较…我最不感到惊讶的结果是：
自定义实现，10个1000x1000matrix乘法的试验：
Took: 15.76542 seconds.BLAS实施，1000×1000matrix乘法的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]; }
我有两个问题：
- 假设一个matrix – matrix乘法：nxm * mxn需要n * n * m次乘法运算，所以在1000 ^ 3或1e9运算的情况下。 BLAS的2.6Ghz处理器如何在1.32秒内完成10 * 1e9操作？ 即使乘法是一个单一的操作，并没有什么别的事情，它应该需要约4秒。
- 为什么我的执行速度如此之慢？
出于很多原因。
首先，Fortran编译器是高度优化的，语言允许它们是这样的。 C和C ++在数组处理方面非常松散（例如指向同一内存区域的指针）。 这意味着编译器不能预先知道该做什么，并且被迫创build通用代码。 在Fortran中，您的案例更加简化，编译器可以更好地控制发生的情况，从而使他优化更多（例如使用寄存器）。
另一件事是Fortran商店专栏，而C存储行数据。 我没有检查你的代码，但要小心你如何执行产品。 在C中，你必须扫描行：这样你扫描你的arrays连续的内存，减lesscaching未命中。 高速caching未命中是低效率的第一个来源。
第三，这取决于你正在使用的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; } } } }
试试吧，我相信你会省下一些东西。
在你的第一个问题上，原因是如果你使用一个简单的algorithm，matrix乘法的比例就是O（n ^ 3）。 有比较好的algorithm。
Robert A. van de Geijn和Enrique S. Quintana-Ortí编着的“编程matrix计算的科学”一书是一个很好的起点。 他们提供免费的下载版本。
BLAS分为三个等级：
-
等级1定义了一组仅对vector进行操作的线性代数函数。 这些function受益于vector化（例如使用SSE）。
-
二级函数是matrix向量运算，例如一些matrix向量乘积。 这些function可以用Level1函数来实现。 但是，如果您可以提供一个专用的实现，使用一些具有共享内存的多处理器体系结构，则可以提高此function的性能。
-
3级function是像matrixmatrix产品一样的操作。 再次，你可以用Level2函数来实现它们。 但是Level3函数对O（N ^ 2）数据执行O（N ^ 3）运算。 因此，如果您的平台具有caching层次结构，那么如果您提供caching优化/caching友好的专用实现，则可以提高性能。 这本书很好地描述。 Level3函数的主要推动来自caching优化。 这个提升大大超过了并行和其他硬件优化的第二个提升。
顺便说一下，大多数（甚至全部）高性能的BLAS实现都不是在Fortran中实现的。 ATLAS在C.中实现。GotoBLAS / OpenBLAS是在C语言中实现的，它在汇编器中是性能关键的部分。 在Fortran中只实现了BLAS的参考实现。 然而，所有这些BLAS实现提供了一个Fortran接口，以便它可以链接到LAPACK（LAPACK从BLAS获得所有的性能）。
优化的编译器在这方面发挥的作用不大（对于GotoBLAS / OpenBLAS编译器根本就不重要）。
恕我直言，没有BLAS实现使用像Coppersmith-Winogradalgorithm或Strassenalgorithm的algorithm。 我不完全确定原因，但这是我的猜测：
- 也许它不可能提供这些algorithm的caching优化实现（即你会失去更多，那么你会赢）
- 这些algorithm在数值上不稳定。 由于BLAS是LAPACK的计算核心，这是一个不可行的方法。
编辑/更新：
关于这个主题的新的和开创性的论文是BLIS论文 。 他们写得非常好。 对于我的讲座“高性能计算的软件基础”，我按照他们的论文实施了matrixmatrix产品。 其实我实现了matrixmatrix产品的几个变种。 最简单的变体完全是用普通的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; }
matrix – matrix乘积的整体性能仅取决于这些回路。 大约99.9％的时间在这里度过。 在其他变体中，我使用内部函数和汇编代码来提高性能。 你可以在这里看到这个教程经历了所有的变种：
ulmBLAS：关于GEMM（matrix – matrix产品）的教程
与BLIS文件一起，了解英特尔MKL等库如何获得这样的性能变得相当容易。 为什么使用行或列主存储无关紧要！
最后的基准在这里（我们称为我们的项目ulmBLAS）：
基准为ulmBLAS，BLIS，MKL，openBLAS和Eigen
另一个编辑/更新：
我还写了一些关于BLAS如何用于求解一个线性方程组的数值线性代数问题的教程：
高性能LU分解
（这个LU分解例如被Matlab用来求解一个线性方程组。
我希望能抽出时间来扩展教程来描述和演示如何在PLASMA中实现LU分解的高度可扩展的并行实现。
好吧，在这里你去编码caching优化并行LU分解
PS：我也做了一些改进uBLAS性能的实验。 这实际上是非常简单的提升（是啊，玩文字:)）的性能的UBLAS：
在uBLAS上的实验 。
这里有一个与BLAZE类似的项目：
BLAZE的实验 。
所以BLAS首先是一个大约50个函数的接口。 接口有许多竞争的实现。
首先我会提到基本上不相关的东西：
- Fortran vs C，没有什么区别
- 先进的matrixalgorithm，如Strassen，实现不使用他们，因为他们没有在实践中帮助
大多数实现或多或less地将每个操作分解为小维matrix或向量操作。 例如，一个大的1000×1000的matrix乘法可能会分解成一个50×50的matrix乘法序列。
这些固定大小的小尺寸操作（称为内核）在CPU特定的汇编代码中使用其目标的多个CPU特性进行硬编码：
- SIMD式指令
- 指令级并行性
- caching意识
而且，这些内核可以使用multithreading（CPU内核）相互并行执行，典型的map-reducedevise模式。
看看最常用的开源BLAS实现ATLAS。 它有许多不同的竞争内核，在ATLAS库构build过程中，它们之间运行竞争（有些甚至是参数化的，所以相同的内核可以有不同的设置）。 它会尝试不同的configuration，然后为特定的目标系统select最佳configuration。
（提示：这就是为什么如果你使用的是ATLAS，那么你最好是专门为你的机器构build和调整库，然后使用预编译的机器。）
首先，matrix乘法的algorithm比你使用的algorithm更高效。
其次，你的CPU一次只能执行多于一条指令。
您的CPU每个周期执行3-4条指令，如果使用SIMD单元，每条指令将处理4个浮点数或2个双精度浮点数。 （当然这个数字也不准确，因为CPU通常每个周期只能处理一个SIMD指令）
第三，你的代码远不是最优的：
- 您正在使用原始指针，这意味着编译器必须假定它们可能是别名。 有编译器专用的关键字或标志，您可以指定告诉编译器，他们没有别名。 或者，您应该使用其他types，而不是原始指针来处理这个问题。
- 你通过对inputmatrix的每一行/列进行一个简单的遍历来颠簸caching。 在移动到下一个块之前，可以使用阻塞在matrix的较小块上执行尽可能多的工作，该matrix适合CPUcaching。
- 对于纯数字任务来说，Fortran几乎是无与伦比的，C ++需要大量的哄骗才能达到类似的速度。 它可以完成，并且有几个库（通常使用expression式模板）来展示它，但这不是微不足道的，它不会发生。
我不清楚BLAS的实现，但matrix乘法有比O（n3）更好的更高效的algorithm。 众所周知的是Strassenalgorithm
第二个问题的大多数论点 – 汇编，分裂成块（但不是N ^ 3以下的algorithm，它们真的过度开发） – 起到了一定的作用。 但是algorithm的低速度主要是由matrix大小和三个嵌套循环的不幸排列造成的。 您的matrix非常大，以至于不能立即放入caching。 您可以重新排列循环，以便尽可能多地在caching中的某一行上完成，这样可以大大减lesscaching刷新（顺便说一下，如果循环遍布块，排列顺序相似，则BTW分为小块有模拟效果。 matrix的模型实现如下。 在我的电脑上，与标准实现相比，它的时间消耗约为1:10（和你的一样）。 换句话说：从来没有按照我们在学校学到的“行时间列”scheme来编程matrix乘法。 在重新安排循环之后，通过展开循环，汇编代码等，可以获得更多的改进。
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_dgemmreplace全部更好（在你的电脑上试试！）。 但更快（1：4）直接调用Fortran库的dgemm_。 我认为这个例程实际上不是Fortran，而是汇编代码（我不知道库中是什么，我没有这个源代码）。 完全不清楚为什么cblas_dgemm没有那么快，因为据我所知，这只是dgemm_的一个包装。
这是一个现实的加速。 有关使用SIMD汇编程序完成C ++代码的示例，请参阅iPhonematrix函数的示例 – 这些函数的速度比C版快8倍，甚至没有“优化”程序集 – 现在还没有pipe道线程是不必要的堆栈操作。
另外你的代码不是“ 限制正确的 ” – 编译器如何知道当它修改C时，它不会修改A和B？
对于MM乘法中的原始代码，大多数操作的内存引用是性能不佳的主要原因。 内存运行速度比caching低100-1000倍。
大部分加速来自MM乘法中的这个三重循环函数的循环优化技术。 使用两种主循环优化技术; 展开和阻止。 关于展开，我们展开最外面的两个循环，并阻止它在caching中的数据重用。 外循环展开有助于在整个操作的不同时间减less对同一数据的内存引用次数，从而暂时优化数据访问。 在特定的数字处阻塞循环索引，有助于将数据保留在caching中。 您可以select优化二级高速caching或三级高速caching。
