寻求最快速的矩阵乘实现

G-Spider

目标： 时间最短是王道! 规则： 1.输入精度： double双精度型（64位），即使用matlab默认的数据精度。 2.输入值： A=rand(n,m),B=rand(m,p) 生成的0~1的两个随机数矩阵(矩阵行列通常在1000位左右)，依然遵循matlab的rand函数 (线性同余伪随机数发生器，可自行实现)，计算C(n,p)=A(n,m)*B(m,p)。 3.时间测试： 从rand开始，到矩阵乘计算结束。 附： 不包括GUI或CUI的显示时间（此显示仅仅确认结果的正确性）。 空间开销不计。指令集不限，程序语言不限，运行平台不限。Good luck!!

gxqcn

目标： ... 3.时间测试： 从rand开始，到矩阵乘计算结束 ... G-Spider 发表于 2010-10-25 23:56

我觉得应是“从rand后开始，到矩阵乘计算结束”，因为rand算法也会有所不同的，虽然影响不大。 我很少涉及浮点运算，也没针对通用矩阵进行优化的经验，所以静观高手表演。

G-Spider

2# gxqcn 郭兄这么说也有道理，对于时间的测试我在写的时候就有所顾虑，矩阵乘核涉及读取数据，计算，返回结果。读取数据的效率会严重影响计算效率，未对齐的数据可能会浪费更多读取时间，我猜想matlab在定义一个大矩阵的时候不会傻瓜似的顺序罗列...我加上rand的时间是希望在数据的预处理上动些手脚，这样可能使矩阵乘核工作的更加有效。我会征求大家意见和建议。 期待高人....

forcal

1# G-Spider 感谢楼主专门出贴讨论这个问题！ 到底是高手，比我在这个帖子21＃的问题2总结的更到位：http://bbs.emath.ac.cn/thread-2709-1-1.html 一直期待国内的高手亮相，老听人说matlab的矩阵算法是最牛的，一直未全信。 做两个测试吧？也不麻烦。 1、从rand后开始，到矩阵乘计算结束。 2、从rand前开始，到矩阵乘计算结束。

forcal

我先献丑，垫个底，助助兴。 1、从rand后开始，到矩阵乘计算结束。 Forcal代码：

1. !using["math","sys"];
2. (:a,b,k,t0)=
3. oo{
4. a=rand[1000,1000], b=rand[1000,1000],
5. t0=clock(),
6. k=a*b //矩阵乘
7. },
8. [clock()-t0]/1000;

复制代码

耗时约2.219秒。 2、从rand前开始，到矩阵乘计算结束。 Forcal代码：

1. !using["math","sys"];
2. (:a,b,k,t0)=
3. t0=clock(),
4. oo{
5. a=rand[1000,1000], b=rand[1000,1000],
6. k=a*b //矩阵乘
7. },
8. [clock()-t0]/1000;

复制代码

耗时约2.328秒。 我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。 只有一个核在运行（多核不会弄）。 可下载OpenFC（绿色免安装，2M）演示以上代码：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

liangbch

可惜我的显卡的GPU不支持double计算，要不我写一个调用CUDA矩阵库的测试程序。

G-Spider

6# liangbch 那就先用单精度型的吧(32位)，实际上我暂时先写的float版本 ，雏形已经有了，使用SSE2系统指令集， 未对B矩阵转置，且A,B的列是4的整数倍的情形下......昨天主要完成了这个模块.... 随机函数使用了c中的rand()改造版本，在种子相同的情形下，得到相同的随机数。 随机参数: 跟踪c版的rand函数： ;....................................................... A(1024,1024)*B(1024,1024) 单精度乘法，A,B种子分别为1和2（单线程在core 2测试已确保结果正确性)，在320ms左右完成，预计加上B的转置和对非4整数倍的处理之后，在500ms内计算完成 直接把部分结果拷到文件中了 当然死算每一个乘积，几乎已到了尽头，虽然用到了矢量.我很不甘心，恐怕这种方式不可能超过matlab的矩阵乘。 刚才看到一篇斯特拉森矩阵乘法的文章我颇有顾虑: 数据从rand a,b到转成0~1之间的浮点，整个过程也只用了15ms左右 接着是计算a*b 乘法的代价过大，只有提供一种高效的算法来减少乘法的运算次数才能有质的飞跃....酝酿中.....

litaoye

矩阵乘法我只知道那个n^log(2,7)的分治方法，难道有类似fft的变换么？

forcal

乘法的代价过大，只有提供一种高效的算法来减少乘法的运算次数才能有质的飞跃....酝酿中..... G-Spider 发表于 2010-10-29 14:24

我在5＃的矩阵乘算法仅使用C/C++实现，如果你也仅使用C/C++实现，能达到或超过我这个速度，那我就帮不上忙了。 否则，使用我这个雕虫小技或许会使速度有更大的提高。

G-Spider

9# forcal 8# litaoye 看来要使用一些高级的算法才能降低时间复杂度. 看一下乘法的计算次数： 普通：0(n^3)......1024^3=1 073 741 824次乘运算 Strassen分治：0(n^2.81)............1024^2.81= 287 701 988次算运算 Don Coppersmith and Shmuel Wmograd改进：0(n^2.367).........1024^2.367= 13 346 887次算运算 ........ 当然在实现上的难度会越来越大............ 最新的矩阵乘算法动态，可参考： Fast and stable matrix multiplication