代码矢量化是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函数计算三重及以上积分。