Wol*_*olf 184 performance matlab cuda matrix-multiplication
我正在使用CUDA,C++,C#和Java进行一些基准测试,并使用MATLAB进行验证和矩阵生成.但是当我乘以MATLAB时,2048x2048甚至更大的矩阵几乎立即成倍增加.
1024x1024 2048x2048 4096x4096
--------- --------- ---------
CUDA C (ms) 43.11 391.05 3407.99
C++ (ms) 6137.10 64369.29 551390.93
C# (ms) 10509.00 300684.00 2527250.00
Java (ms) 9149.90 92562.28 838357.94
MATLAB (ms) 75.01 423.10 3133.90
只有CUDA具有竞争力,但我认为至少C++会有点接近并且不会60x慢.
所以我的问题是 - MATLAB如何快速地完成它？
C++代码:
float temp = 0;
timer.start();
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] * matice2[m][k];
}
matice3[j][k] = temp;
}
}
timer.stop();
编辑:我也不知道如何考虑C#结果.该算法是一样的为C++和Java,但有一个巨大的跳跃2048的1024？
Edit2:更新了MATLAB和4096x4096结果
rev*_*eer 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算法是错误的,因为上面提到的硬件考虑因素,这两种算法在实践中都不可实现.
Edr*_*ric 83
这是我在使用特斯拉C2070的机器上使用MATLAB R2011a + Parallel Computing Toolbox的结果:
>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.
MATLAB使用高度优化的库进行矩阵乘法,这就是普通MATLAB矩阵乘法如此之快的原因.该gpuArray版本使用MAGMA.
更新使用R2014a与特斯拉K20C的机器上,而新的timeit和gputimeit功能:
>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
0.0324
>> gputimeit(@()gA*gA)
ans =
0.0022
在具有16个物理内核和Tesla V100的WIN64机器上使用R2018b进行更新:
>> timeit(@()A*A)
ans =
0.0229
>> gputimeit(@()gA*gA)
ans =
4.8019e-04
Jon*_*ore 19
Matlab在不久前收录了LAPACK,所以我假设他们的矩阵乘法使用至少那么快的东西.LAPACK源代码和文档随时可用.
你也可以看看转到和van de Geijn的论文"高性能矩阵乘法的解剖"在http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf
小智 9
在进行矩阵乘法时,使用花费时间的朴素乘法O(n^3).
存在矩阵乘法算法O(n^2.4).这意味着n=2000您的算法需要的计算量是最佳算法的100倍.
您应该检查维基百科页面中的矩阵乘法,以获取有关实现它的有效方法的更多信息.
根据您的Matlab版本,我相信它可能已经在使用您的GPU了.
另一件事; Matlab会跟踪矩阵的许多属性; 它的对角线,hermetian等等,并专门研究其基于此的算法.也许它的专业化基于您传递的零矩阵,或类似的东西？也许它正在缓存重复的函数调用,这会扰乱你的时间？也许它优化了重复使用的矩阵产品？
为了防止发生这种情况,请使用随机数字矩阵,并确保通过将结果打印到屏幕或磁盘或其他某些部分来强制执行.