极限测试之Matlab与Forcal代码矢量化

代码矢量化是matlab的特色，但这点似乎不难实现。代码矢量化的优势并不明显，通过一个例子说明。

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

Matlab代码：

tic [x,y]=meshgrid(0:0.0011:1,1:0.0009:2); s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))) toc s = 1.0086e+006 Elapsed time is 0.561108 seconds.

Forcal代码：

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

结果：
1008606.64947441
0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。

Matlab代码：

tic s=0; x = 0; y = 1; while x<1 while y<2; s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); y = y+0.0009; end x = x+0.0011; y = 1; end s toc s = 1.0086e+006 Elapsed time is 0.933513 seconds.

Forcal代码：

mvar: t=sys::clock(); s=0,x=0, while{x<=1, //while循环算法； y=1, while{y<=2, s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), y=y+0.0009 }, x=x+0.0011 }, s; [sys::clock()-t]/1000;

结果：
1008606.64947441
0.734 //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

Forcal中还有一个函数sum专门进行这种计算：

mvar: t=sys::clock(); f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); sum["f",0,1,0.0011 : 1,2,0.0009]; [sys::clock()-t]/1000;

结果：
1008606.64947441
0.719 //时间，秒

＝＝＝＝＝＝＝＝＝＝＝＝＝

＝＝＝＝＝＝＝＝＝＝＝＝＝

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。

matlab代码：

f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); tic [x,y]=meshgrid(0:0.0011:1,1:0.0009:2); sum(sum(arrayfun(f,x,y))) toc ans = 1.0086e+006 Elapsed time is 7.339923 seconds.

Forcal代码：

!using["math","sys"]; mvar: f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); t=clock(), oo{ ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y], arrayfun[HFor("f"),x,y].Sum[] }; [clock()-t]/1000;

结果：
1008606.64947441
0.735秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。