寻找最高效的矢量化算法

forcal

10# G-Spider
我的理解是，矢量化就是对线性排列的数据（数组）整体操作，这样可以提高速度。
matlab常常将一些近乎零空间的操作展开：将操作数分布到数组中以实现所谓的矢量化操作，这就是所谓的以空间换时间。

举例来说：

1. //用C++代码描述为：
   
2. s=0.0;  
   
3. for(x=0.0;x<=1.0;x=x+0.0011)  
   
4. {
   
5.   for(y=1.0;y<=2.0;y=y+0.0009)
   
6.   {
   
7.   s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
8.   }
   
9. }  

复制代码

matlab代码：

1. %file speedtest.m
   
2. function speedtest
   
3. format long
   
4. tic
   
5. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   
6. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   
7. toc

复制代码

大家注意一下meshgrid函数，该函数以两个数组[0:0.0011:1]和[1:0.0009:2]作为输入参数（数组大小分别是910和1112），输出两个大小为910×1112的数组x和y。这就是所谓的代码展开。

然后，sum函数中y/2将数组y所有元素除以2，得到一个同样大小的新数组参与后续计算；x.^(y/2)将数组x与数组(y/2)的对应元素做乘方运算，得到一个同样大小的新数组参与后续计算；... ...，以此类推。最后得到的数组（此例中，该数组实际是一个矩阵）通过两个sum求和。内层的sum先对列求和，外层的sum再对行求和。最后得到一个数。

通过这样处理，运行速度和上面的C/C++代码相当。然而，matlab对这里面涉及的线性运算又做了优化，使用了C/C++代码以外的手段，故效率超过了C/C++。

我这个帖子的用意就是让大家讨论这个C/C++代码以外的手段，即1＃我的问题。

顺便给出此例的Forcal代码：

1. !using["math","sys"];
   
2. mvar:
   
3. t=clock(),
   
4. oo{
   
5.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
   
6.   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn(2)))))),0]
   
7. };
   
8. [clock()-t]/1000;

复制代码

Forcal代码仅达到了C/C++代码的速度（差别很小）。

大家集思广益，看能否将这个问题解决了。欢迎任何推荐的资料和建议，更希望能看到得以实现的代码。

forcal

还有一个矩阵乘的例子，参考：http://forum.simwe.com/archiver/tid-953713.html

lin2009 给出的测试结果似乎说明Forcal的矩阵乘与matlab效率差不多，但bainhome 的测试结果似乎又说明有几倍的差距，我看不懂这个，希望有人给分析一下？

G-Spider

可以说，你跟我所看到的是一个不同层次的加速。"当高层次的理论趋于成熟的时候，我们开始抱怨硬件。"
比如，对于你说的"sum函数中y/2将数组y所有元素除以2"，你可能忽略了对这个计算的加速，
而对于SIMD正是想对这一方面的底层加速，NVIDIA的GPU平台也是对这方面的考虑，
即将诞生的AVX指令集也是如此，这个是不容忽视的。我想这也不是使用C/C++代码的优化手段。
为什么常听人说要进行人工优化？或许很多人会说，我只要加个优化编译选项就好了，...这是远远不够的。
比如对sum函数中y/2将数组y所有元素除以2，是针对一个数一个数进行除以2（对于.x86常用指令只能做到这个点）还是同时对几个数同时除以2（SSE指令集可以做到），这种细节上的加速当然会影响提速。
我想将来关于计算库（当前也有很多这样的尝试），通常会有针对不同指令集的核心库优化，从.X86到SSE2到SS3或更高版本到GPU。
对这些指令集的编程会同时实现，当侦测到支持GPU，用调用GPU算核版本，否则看是否支持SSE3，如果支持就调用SSE3实现的算核版本，依次....
刚看到有人说：“matlab在不同的计算机上运行效率差别怎么这么大？”我想至少有这方面的考虑....
这方面有个经典的例子：FFFF - Fast Floating Fractal Fun  by pacca

就矩阵乘，我当时也尝试用SSE实现了，由于没有进一步的需求，所以也没深入下去。

gxqcn

看来 G-Spider 是同道之人，深谙指令集加速技术。
通过自动侦测CPU支持的指令集，然后自动切换成最高效的指令，从而自动调整算法以充分发挥CPU之功效。
我现在也在等AVX指令集的普及。

BTW，据说64位的编译器不支持内嵌汇编，不知道是否是真的？

无心人

郭兄，何必要求内嵌汇编？
为啥不写独立汇编文件呢？
===================
另外，利用C++的模板，有人研究出了延迟计算技术
利用模板，将要计算的矩阵表达式，展开成为最终式子
而不是一步一算

这就是矩阵运算速度快的原因

gxqcn

我以前还没有写独立汇编文件的经验。

模板，是不是只在编译期起作用啊？在运行期无法自动展开吧？

gxqcn

我觉得楼主这个问题需要从如下几点考虑：
1、精度需求：float、double 虽都为浮点，但精度不同，所需的代价肯定不同，你需求的是什么？MatLab 的精度是怎样的？
2、是否可循环展开？因为现在CPU可以多指令并行（前后无关联）
3、是否可用先进的指令集？
4、是否可多核并行？

forcal

可以说，你跟我所看到的是一个不同层次的加速。"当高层次的理论趋于成熟的时候，我们开始抱怨硬件。"
比如，对于你说的"sum函数中y/2将数组y所有元素除以2"，你可能忽略了对这个计算的加速，
而对于SIMD正是想对这 ...
G-Spider 发表于 2010-10-24 18:42

呵呵，我是按我对matlab以空间换时间的代码矢量化的理解来解释我在1＃提出的问题，不但没有忽略“sum函数中y/2将数组y所有元素除以2”这个计算的加速，而且就是要得到这个加速。

虽然这样，以空间换时间的做法是否可取，仍然是一个有争议的问题，不过我们先不谈这个。

“我想将来关于计算库（当前也有很多这样的尝试），通常会有针对不同指令集的核心库优化，从.X86到SSE2到SS3或更高版本到GPU。....
就矩阵乘，我当时也尝试用SSE实现了，由于没有进一步的需求，所以也没深入下去。”

很期待这样的尝试哦，若实现了这样的算法库，可以专门为Forcal这个不起眼的脚本设计一个“Forcal版G-Spider 计算库”，很期待！呵呵。

forcal

利用C++的模板，有人研究出了延迟计算技术
利用模板，将要计算的矩阵表达式，展开成为最终式子
而不是一步一算

这就是矩阵运算速度快的原因 无心人 发表于 2010-10-25 08:49

这种“延迟计算技术”能否详细说一下？我想很多像我这样的外行愿闻其详！谢谢！

forcal

我觉得楼主这个问题需要从如下几点考虑：
1、精度需求：float、double 虽都为浮点，但精度不同，所需的代价肯定不同，你需求的是什么？MatLab 的精度是怎样的？
2、是否可循环展开？因为现在CPU可以多指令并行（前后无关联）
3、是否可用先进的指令集？
4、是否可多核并行？ gxqcn 发表于 2010-10-25 09:48

感谢gxqcn 的意见！
1、一般，MatLab 使用double 计算。
2和3我不懂，呵呵。
4、我在1＃的问题应该可以多核并行吧？