一个很好的起点是Robert A.货车de Geijn和Enrique S.昆塔纳-奥蒂他们提供免费下载版本。
BLAS分为三个级别：

• 级别1定义了一组仅对向量进行操作的线性代数函数。这些函数受益于向量化（例如，使用SIMD，如SSE）。
• 2级函数是矩阵向量运算，例如一些矩阵向量积。这些功能可以按照1级功能来实现。但是，如果您能够提供一个专用的实现，利用一些具有共享内存的多处理器体系结构，那么您就可以提高这些函数的性能。
• 3级函数是类似矩阵-矩阵乘积的运算。同样，您可以根据第2级函数来实现它们。但是3级函数在O（N^2）数据上执行O（N^3）操作。因此，如果您的平台具有缓存层次结构，那么如果您提供 * 缓存优化/缓存友好 * 的专用实现，则可以提高性能。这在书中有很好的描述。Level 3函数的主要提升来自缓存优化。这种提升明显超过了并行和其他硬件优化的第二次提升。

顺便说一下，大多数（甚至所有）高性能BLAS实现都不是用Fortran实现的。ATLAS是用C语言实现的。GotoBLAS/OpenBLAS是用C实现的，其性能关键部分是用汇编实现的。在Fortran中只实现了BLAS的参考实现。但是，所有这些BLAS实现都提供了Fortran接口，因此可以与LAPACK链接（LAPACK从BLAS获得所有性能）。
优化的编译器在这方面起着次要的作用（对于GotoBLAS/OpenBLAS，编译器根本不重要）。
恕我直言，无BLAS实现使用像Coppersmith-Winograd算法或斯特拉森算法的算法。可能的原因是：

• 也许不可能提供这些算法的缓存优化实现（即你输的比赢的多）
• 这些算法在数值上是不稳定的。由于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;
    }

矩阵-矩阵乘积 only 的整体性能取决于这些循环。99.9%的时间都在这里度过。在其他变体中，我使用了intrinsic和汇编代码来提高性能。你可以在这里看到教程中的所有变体：
ulmBLAS: Tutorial on GEMM (Matrix-Matrix Product)
结合BLIS论文，我们可以很容易地理解像英特尔MKL这样的库如何获得这样的性能。以及为什么使用行主存储还是列主存储并不重要！
最后的基准测试在这里（我们称我们的项目为ulmBLAS）：
ulmBLAS、BLIS、MKL、openBLAS和Eigen的基准测试
编辑/更新：
我还写了一些关于如何将BLAS用于数值线性代数问题的教程，例如求解线性方程组：
High Performance LU Factorization
(This LU因式分解例如由Matlab用于求解线性方程组。
我希望能抽出时间来扩展本教程，以描述和演示如何实现PLASMA中LU分解的高度可扩展并行实现。
好的，给你：Coding a Cache Optimized Parallel LU Factorization
P.S.：我也做了一些关于提高uBLAS性能的实验。实际上，提高uBLAS的性能非常简单（是的，玩文字游戏：）：
Experiments on uBLAS。
这里有一个类似的项目与BLAZE：
Experiments on BLAZE。

首先，BLAS只是一个包含大约50个函数的接口。该接口有许多相互竞争的实现。
首先，我会提到一些基本上不相关的事情：

• Fortran vs C，没有区别
• 高级矩阵算法，如斯特拉森，实现不使用它们，因为它们在实践中没有帮助

大多数实现以或多或少明显的方式将每个操作分解为小维度矩阵或向量操作。例如，大的1000x1000矩阵乘法可以分解成50x50矩阵乘法的序列。
这些固定大小的小维度操作（称为内核）使用其目标的几个CPU功能在CPU特定的汇编代码中硬编码：

• SIMD风格的指令
• 指令级并行
• 缓存感知

此外，在典型的Map缩减设计模式中，这些内核可以使用多个线程（CPU内核）相对于彼此并行执行。
看看ATLAS，它是最常用的开源BLAS实现。它有许多不同的竞争内核，在ATLAS库构建过程中，它在它们之间运行竞争（有些甚至是参数化的，因此同一内核可以有不同的设置）。它会尝试不同的配置，然后为特定的目标系统选择最佳配置。
(Tip：这就是为什么如果您使用ATLAS，那么最好为您的特定机器手动构建和调优库，然后使用预构建的库。

首先，有比您正在使用的算法更有效的矩阵乘法算法。
第二，CPU一次可以执行多条指令。
CPU每个周期执行3-4条指令，如果使用SIMD单元，每条指令处理4个浮点数或2个双精度数。（当然这个数字也不准确，因为CPU通常每个周期只能处理一条SIMD指令）
第三，你的代码远远不是最佳的：

• 您使用的是原始指针，这意味着编译器必须假设它们可能有别名。您可以指定特定于编译器的关键字或标志来告诉编译器它们不使用别名。或者，您应该使用原始指针以外的其他类型，这样可以解决这个问题。
• 通过对输入矩阵的每一行/每一列执行简单的遍历，您会该高速缓存变得混乱不堪。您可以使用分块在矩阵的较小块上执行尽可能多的工作，该块适合CPU缓存，然后再移动到下一个块。
• 对于纯数值任务，Fortran几乎是无与伦比的，而C++需要大量的哄骗才能达到类似的速度。这是可以做到的，并且有一些库演示了它（通常使用表达式模板），但它并不简单，而且它不会“突然”发生。

我不知 prop 体的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_dgemm没有那么快，因为据我所知，它只是dgemm_的 Package 器。

这是一个现实的加速。关于SIMD汇编器在C++代码上可以做什么的示例，请参见一些示例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++在数组处理方面非常宽松（例如指针指向同一存储器区域的情况）。这意味着编译器无法提前知道要做什么，并且被迫创建泛型代码。在Fortran中，你的用例更加流线化，编译器可以更好地控制发生的事情，允许他优化更多（例如：使用寄存器）。
另一件事是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;
           }
     }
  }
}

试试看，我相信你会保存一些东西。
对于你的第一个问题，原因是如果你使用一个平凡的算法，矩阵乘法的规模为O（n^3）。有一些算法是scale much better。