reverse_engi..

163

这种问题反复出现,应该比Stackoverflow上的"Matlab使用高度优化的库"或"Matlab使用MKL"一次更清晰地回答.

历史:

矩阵乘法(与矩阵向量,向量 - 向量乘法和许多矩阵分解一起)是线性algrebra中最重要的问题.从早期开始,工程师就一直在用计算机解决这些问题.

我不是历史专家,但显然那时候,每个人都只用简单的循环重写了他的Fortran版本.然后出现了一些标准化,识别出需要解决的大多数线性代数问题的"内核"(基本例程).然后将这些基本操作标准化为:基本线性代数子程序(BLAS).然后,工程师可以在他们的代码中调用这些经过良好测试的标准BLAS例程,使他们的工作变得更加容易.

BLAS:

BLAS从1级(定义标量矢量和矢量矢量运算的第一个版本)演变为2级(矢量矩阵运算)到3级(矩阵矩阵运算),并提供越来越多的"内核",使标准化更多以及更多基本的线性代数运算.最初的Fortran 77实现仍可在Netlib的网站上获得.

为了更好的表现:

因此,多年来(特别是在BLAS 1级和2级版本之间:80年代早期),随着向量操作和缓存层次结构的出现,硬件发生了变化.这些演进使得有可能大大提高BLAS子程序的性能.然后不同的供应商出现了BLAS例程的实现,这些例程越来越高效.

我不知道所有的历史实现(当时我还没出生或是个孩子),但是最着名的两个是在21世纪初出现的:英特尔MKL和GotoBLAS.您的Matlab使用的是英特尔MKL,这是一款非常优秀的优化BLAS,它解释了您所看到的出色性能.

Matrix乘法的技术细节:

那么为什么Matlab(MKL)如此快dgemm(双精度一般矩阵 - 矩阵乘法)呢？简单来说:因为它使用矢量化和良好的数据缓存.更复杂的术语:请参阅Jonathan Moore提供的文章.

基本上,当您在所提供的C++代码中执行乘法运算时,您根本不需要缓存.因为我怀疑你创建了一个指向行数组的指针数组,所以你在内部循环中访问"matice2"的第k列:matice2[m][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的大小(转换后备缓冲区,长话短说:可以有效缓存的内容)定义的矩阵的块执行乘法,以便它们流式传输到处理器确切地说它可以处理的数据量.另一方面是矢量化,它们使用处理器的矢量化指令来获得最佳指令吞吐量,而您无法通过跨平台C++代码实现这一点.

最后,人们声称这是因为Strassen或Coppersmith-Winograd算法是错误的,因为上面提到的硬件考虑因素,这两种算法在实践中都不可实现.

我刚刚观看了Scott Meyers的视频，介绍了缓存大小的重要性以及使数据适合缓存行大小的问题，以及在源中没有共享数据但最终在硬件上共享数据的多线程解决方案可能会遇到的问题/ core-thread level：https://youtu.be/WDIkqP4JbkE (2认同)