matlab指针矩阵乘法,为什么MATLAB在矩阵乘法方面如此之快？

这种问题反复出现，应该比“MATLAB使用高度优化的库”或“MATLAB使用MKL”一次堆栈溢出更清楚地回答。

历史记录：

矩阵乘法(连同矩阵向量、向量乘法和许多矩阵分解)是线性代数中最重要的问题。工程师们从早期就开始用计算机解决这些问题。

我不是历史方面的专家，但显然在那时，每个人都只是用简单的循环重写他的FORTRAN版本。随后出现了一些标准化，并确定了大多数线性代数问题需要解决的“核心”(基本例程)。这些基本操作随后在一个称为：基本线性代数子程序(BLAS)的规范中标准化。工程师可以在他们的代码中调用这些标准的、经过良好测试的BLAS例程，使他们的工作更加容易。

布拉斯：

BLAS从1级(定义标量向量和向量向量运算的第一个版本)发展到2级(向量矩阵运算)到3级(矩阵矩阵运算)，并提供了越来越多的“核”，使越来越多的基本线性代数运算标准化。最初的FORTRAN 77实现仍然可以在Netlib's website上使用。

提高性能：

因此，多年来(特别是在BLAS 1级和2级版本之间：80年代早期)，随着矢量操作和缓存层次结构的出现，硬件发生了变化。这些演进使得大幅度提高BLAS子例程的性能成为可能。随后，不同的供应商推出了他们的BLAS例程的实现，这些例程的效率越来越高。

我不知道所有的历史实现(那时我还没有出生或是个孩子)，但最著名的两个实现出现在21世纪初：英特尔MKL和GotoBLAS。您的Matlab使用Intel MKL，这是一个非常好的、优化的BLAS，这解释了您看到的出色性能。

矩阵乘法的技术细节：

那么，为什么Matlab(MKL)在dgemm(双精度通用矩阵乘法)上速度如此之快呢？简单来说：因为它使用矢量化和良好的数据缓存。用更复杂的术语来说：请参阅Jonathan Moore提供的article。

基本上，当你在你提供的C++代码中执行乘法时，你根本就不喜欢缓存。因为我怀疑您创建了一个指向行数组的指针数组，所以您在内部循环中对“matice2”：matice2[m][k]的第k列的访问非常慢。实际上，当您访问matice2[0][k]时，必须获得矩阵数组0的第k个元素。然后在下一次迭代中，必须访问matice2[1][k]，这是另一个数组(数组1)的第k个元素。然后在下一次迭代中访问另一个数组，依此类推。。。由于整个矩阵matice2不能放在最高的缓存中(它的8*1024*1024字节大)，程序必须从主内存中获取所需的元素，从而浪费大量时间。

如果您只是转置矩阵，以便访问位于连续的内存地址中，那么您的代码将运行得更快，因为现在编译器可以同时将整行加载到缓存中。请尝试此修改版本：timer.start();

float temp = 0;

//transpose matice2

for (int p = 0; p < rozmer; p++)

{

for (int q = 0; q < rozmer; q++)

{

tempmat[p][q] = matice2[q][p];

}

}

for(int j = 0; j < rozmer; j++)

{

for (int k = 0; k < rozmer; k++)

{

temp = 0;

for (int m = 0; m < rozmer; m++)

{

temp = temp + matice1[j][m] * tempmat[k][m];

}

matice3[j][k] = temp;

}

}

timer.stop();

所以您可以看到缓存局部性是如何大大提高代码的性能的。现在，真正的dgemm实现将这一点扩展到了一个非常广泛的层次：它们在由TLB(Translation lookaside buffer，long story short:what can effective be cached)大小定义的矩阵块上执行乘法，以便它们将它能够处理的数据量准确地流到处理器。另一个方面是矢量化，它们使用处理器的矢量化指令来获得最佳的指令UpPoT，你不能从你的跨平台C++代码中真正做到。

最后，人们声称这是因为斯特拉森或铜匠-温诺格拉德算法是错误的，这两种算法在实践中都无法实现，因为硬件考虑上述。