出于好奇,我决定对我自己的矩阵乘法函数与BLAS实现进行基准测试。我不得不说,我对这个结果并不感到惊讶:
自定义实现,10次1000 x1000矩阵乘法试验:
Took: 15.76542 seconds.
BLAS实施,10次1000 x1000矩阵乘法试验:
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];
}
我有两个问题:
1.假设一个矩阵-矩阵乘法说:nxm * mxn需要nnm次乘法,因此在1000^3或1 e9次运算以上的情况下。在我的2.6Ghz处理器上,BLAS如何在1.32秒内完成10* 1 e9操作?即使乘法是一个单一的操作,没有其他事情要做,它应该需要~4秒。
1.为什么我的实现速度这么慢?
8条答案
按热度按时间ljo96ir51#
一个很好的起点是Robert A.货车de Geijn和Enrique S.昆塔纳-奥蒂他们提供免费下载版本。
BLAS分为三个级别:
顺便说一下,大多数(甚至所有)高性能BLAS实现都不是用Fortran实现的。ATLAS是用C语言实现的。GotoBLAS/OpenBLAS是用C实现的,其性能关键部分是用汇编实现的。在Fortran中只实现了BLAS的参考实现。但是,所有这些BLAS实现都提供了Fortran接口,因此可以与LAPACK链接(LAPACK从BLAS获得所有性能)。
优化的编译器在这方面起着次要的作用(对于GotoBLAS/OpenBLAS,编译器根本不重要)。
恕我直言,无BLAS实现使用像Coppersmith-Winograd算法或斯特拉森算法的算法。可能的原因是:
编辑/更新:
这个主题的新的和开创性的论文是BLIS papers。它们写得特别好。在我的演讲“高性能计算的软件基础”中,我实现了他们论文之后的矩阵-矩阵乘积。实际上,我实现了矩阵-矩阵乘积的几个变体。最简单的变体完全用纯C语言编写,只有不到450行代码。所有其他变体都只是优化了循环
矩阵-矩阵乘积 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。
jmo0nnb32#
首先,BLAS只是一个包含大约50个函数的接口。该接口有许多相互竞争的实现。
首先,我会提到一些基本上不相关的事情:
大多数实现以或多或少明显的方式将每个操作分解为小维度矩阵或向量操作。例如,大的1000x1000矩阵乘法可以分解成50x50矩阵乘法的序列。
这些固定大小的小维度操作(称为内核)使用其目标的几个CPU功能在CPU特定的汇编代码中硬编码:
此外,在典型的Map缩减设计模式中,这些内核可以使用多个线程(CPU内核)相对于彼此并行执行。
看看ATLAS,它是最常用的开源BLAS实现。它有许多不同的竞争内核,在ATLAS库构建过程中,它在它们之间运行竞争(有些甚至是参数化的,因此同一内核可以有不同的设置)。它会尝试不同的配置,然后为特定的目标系统选择最佳配置。
(Tip:这就是为什么如果您使用ATLAS,那么最好为您的特定机器手动构建和调优库,然后使用预构建的库。
fafcakar3#
首先,有比您正在使用的算法更有效的矩阵乘法算法。
第二,CPU一次可以执行多条指令。
CPU每个周期执行3-4条指令,如果使用SIMD单元,每条指令处理4个浮点数或2个双精度数。(当然这个数字也不准确,因为CPU通常每个周期只能处理一条SIMD指令)
第三,你的代码远远不是最佳的:
pgx2nnw84#
我不知 prop 体的BLAS实现,但有更有效的算法矩阵乘法,具有优于O(n3)的复杂度。一个众所周知的例子是Strassen Algorithm
5hcedyr05#
第二个问题的大多数论点--汇编程序,拆分成块等。(但不是少于N^3的算法,它们确实过度开发了)--发挥作用。但是你的算法的低速度本质上是由矩阵大小和三个嵌套循环的不幸安排造成的。您的矩阵太大了,以至于它们无法立即放入缓存内存。你可以重新安排循环,这样尽可能多的将在缓存中的一行上完成,这种方式大大减少了缓存刷新(顺便说一句,分割成小块具有类似的效果,最好是在块上的循环以类似的方式排列)。以下是方阵的模型实现。在我的计算机上,与标准实现(如您的)相比,其时间消耗约为1:10。换句话说:永远不要沿着我们在学校里学过的“行乘列”的方式来编程矩阵乘法。在重新安排循环之后,通过展开循环、汇编代码等获得更多的改进。
还有一句话在我的计算机上,这个实现甚至比用BLAS例程cblas_dgemm替换所有的更好(在您的计算机上试试!)。但直接调用Fortran库的dgemm_要快得多(1:4)。我认为这个例程实际上不是Fortran,而是汇编代码(我不知道库中有什么,我没有源代码)。我完全不清楚为什么cblas_dgemm没有那么快,因为据我所知,它只是dgemm_的 Package 器。
yduiuuwa6#
这是一个现实的加速。关于SIMD汇编器在C++代码上可以做什么的示例,请参见一些示例iPhone matrix functions-这些比C版本快8倍以上,甚至不是“优化”的汇编-还没有管道衬里,并且有不必要的堆栈操作。
你的代码也不是“restrict correct”--编译器怎么知道当它修改C时,它没有修改A和B?
rsaldnfx7#
对于MM乘法的原始代码,大多数操作的内存引用是导致性能差的主要原因。内存的运行速度比缓存慢100-1000倍。
大部分的加速来自于对MM乘法中的三重循环函数采用循环优化技术。使用了两种主要的循环优化技术;展开和阻挡。关于展开,我们展开最外面的两个循环,并将其阻塞,以便在缓存中重用数据。外循环展开通过减少在整个操作期间在不同时间对相同数据的存储器引用的数量来帮助暂时地优化数据访问。在特定的数字处阻塞循环索引,有助于将数据保留在缓存中。您可以选择优化二级缓存或三级缓存。
https://en.wikipedia.org/wiki/Loop_nest_optimization
jk9hmnmh8#
有很多原因。
首先,Fortran编译器是高度优化的,语言允许它们这样做。C和C++在数组处理方面非常宽松(例如指针指向同一存储器区域的情况)。这意味着编译器无法提前知道要做什么,并且被迫创建泛型代码。在Fortran中,你的用例更加流线化,编译器可以更好地控制发生的事情,允许他优化更多(例如:使用寄存器)。
另一件事是Fortran按列存储数据,而C按行存储数据。我还没有检查你的代码,但要小心你如何执行产品。在C中,你必须逐行扫描:通过这种方式,您可以沿着连续内存扫描阵列,从而减少该高速缓存未命中。高速缓存未命中是效率低下的第一个原因。
第三,这取决于您使用的blas实现。有些实现可能是用汇编语言编写的,并针对您正在使用的特定处理器进行了优化。netlib版本是用fortran 77编写的。
此外,您正在执行大量操作,其中大多数是重复和冗余的。所有那些获得索引的乘法对于性能是有害的。我真的不知道在BLAS中是如何做到这一点的,但有很多技巧可以防止昂贵的操作。
例如,您可以用这种方式重新编写代码
试试看,我相信你会保存一些东西。
对于你的第一个问题,原因是如果你使用一个平凡的算法,矩阵乘法的规模为O(n^3)。有一些算法是scale much better。