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

forcal

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

1. //用C++代码描述为：
   , d3 d  F4 O# [0 U4 y
2. s=0.0; + ^& {/ q1 `' ]
   
3. for(x=0.0;x<=1.0;x=x+0.0011) ( s8 ~/ j0 z# F( C0 e5 P
   
4. {
   2 Z$ ?/ L3 P1 k* F
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   ( d: O\" p# T0 [- s' f4 |( _
6.    {
   # G) o# p4 q- }/ ]* `( v8 E
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));9 Z, `1 @9 i# r; g. N
   
8.    }
   5 G; K* j* w' Y
9. }

复制代码

Matlab代码：

1. tic6 G6 d7 w) [5 S& ~( q( \
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   5 _4 b+ \& N- y1 N8 o$ q- L7 n
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ! z7 C' |5 v; I  M
4. toc
   & p0 p5 f  s% d+ K
5. 
   ( l\" ~7 g' e9 Q6 u\" \
6. s =
   0 Z5 G* M2 w2 }/ \1 u) @1 G
7. 
   4 |; }, z0 l; U' X. H! u, U
8.   1.0086e+006
   # b, m. f5 E; h$ |5 i5 F
9. ; [7 x1 Q% D  \+ n2 v2 I, o
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];. Z! _. @- O1 k- |: I
   
4. mvar:
   
5. , W$ b. r4 \& z9 M+ ?
6. t=clock(),. U3 F  Q9 b# P/ @  F/ i
   
8. oo{
   
9. 4 e7 {$ V# A0 Y7 g9 }0 Y+ O# ^
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: h% S2 |& {1 O- I/ C$ G
    
12.     Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
    
13. . Q; P5 F2 e' z0 z% o
14. };
    
15. % M% a% e1 K$ Q7 ?& b; i' i0 h& t
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">. Z! _. @- O1 k- |: I</font><br /> mvar:<br /> <span style=\\"display:none\\">, W$ b. r4 \& z9 M+ ?</span> t=clock(),<font class=\\"jammer\\">. U3 F  Q9 b# P/ @  F/ i</font><br /> oo{<br /> <span style=\\"display:none\\">4 e7 {$ V# A0 Y7 g9 }0 Y+ O# ^</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: h% S2 |& {1 O- I/ C$ G</font><br />     Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]<br /> <span style=\\"display:none\\">. Q; P5 F2 e' z0 z% o</span> };<br /> <span style=\\"display:none\\">% M% a% e1 K$ Q7 ?& b; i' i0 h& t</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641

Forcal比Matlab稍慢些。
! a/ p, F4 J/ h' T! w0 U! a
; _+ y5 Y* z3 H9 S－－－－－－－－－－
; n  z' }# A# p/ s5 T) [2 ?! Z- A
5 d: `1 Z7 p' Y. z再看循环效率。
& h% P( g" U* Q4 B+ t
, I; ~0 G" [% i, F4 X6 o6 QMatlab代码：

1. tic
   2 {& j9 L  ?6 c# H5 K
2. s=0;3 o) h- I/ y; ~! P
   
3. x = 0;1 ?+ E% n, a\" y; t
   
4. y = 1;
   ; ^; x8 H# T9 t3 A. u2 [1 L
5. while x<1    + Y' m. Y. Z/ F. X# n
   
6.     while y<2;        ) d2 c+ n$ N\" p
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" r: L. ~# Q0 p' j
8.         y = y+0.0009;        
   . I4 `' G) T/ V( @
9.     end
   ' J0 y! I  s2 R6 F7 L3 S\" z3 S
10.     x = x+0.0011;, x& k- \$ m5 ]
    
11.     y = 1;
    7 \# U. p4 ~  S\" g' G/ l) j
12. end
    6 ^  r1 k* K: Y0 a! U
13. s9 Q) {$ o# Y: |  l
    
14. toc) B* S6 R\" q& N2 W7 \
    
15. 8 {6 R$ h* ~3 r
    
16. s =
    6 Y; V$ [) j# c+ y5 D
17. 
    ) B# Y\" X! j( \\" N! D
18.   1.0086e+006, l5 j- V( M% N) V
    
19. 4 m$ G. P! J5 u' E# @3 m5 U& w
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:2 h: Y: E! y7 a6 z, c: T
   
2. t=sys::clock();
   ! p% @\" F! `+ j! L\" h9 Y6 T
3. s=0,x=0, # w& e( t3 g' K: I# `- o1 o! U8 q5 r: t
   
4. while{x<=1, //while循环算法； 9 @) ~, o7 E8 I9 _2 C
   
5.     y=1, : A: U3 t  ~/ S. V& y& d4 J
   
6.     while{y<=2, ! L. ^3 X! D7 X( l1 \5 b( [
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 3 u\" u7 c; {+ X7 j; K: D
   
8.         y=y+0.0009
   & o$ I; ]8 E+ J& V& r) R1 Z: c) s, t
9.     }, ) n+ u7 v  N6 s
   
10.     x=x+0.0011
    5 U1 @5 j9 M\" w* Z* g
11. }, # K3 x\" \) R5 n
    
12. s;# H& C9 P$ t1 F/ O6 K
    
13. [sys::clock()-t]/1000;

复制代码

结果：
! G! G; }' D, q* ^/ o1008606.64947441
6 N8 ~5 U. V( {0.734   //时间，秒
" R, J0 D- R* S& p$ ^! N
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

" P- D$ p% Y& I8 l6 `+ M/ y$ j5 x5 f9 w－－－－－－－

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

1. mvar:; }4 `; J+ j6 c4 V+ H4 u
   
2. t=sys::clock();
   : t; u& T, v9 V: H
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . E: S0 w5 p8 X! L6 M
4. sum["f",0,1,0.0011 : 1,2,0.0009];: p) `5 c+ r* I$ {
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
2 |" ^' O0 d: ?: a2 U) q
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 \. ~' H5 s  [2 b, \- f
2. tic
   % k0 V! F/ j7 @9 S% P
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);; n+ x$ l2 N/ \3 V  k0 C8 i- c1 D- I
   
4. sum(sum(arrayfun(f,x,y)))
   ( i+ a4 V* P2 E: W; p. Y: _- S
5. toc
   2 i, W6 ?4 K' O7 S' s' N
6. 0 ~# J& X) i) f8 q0 E
   
7. ans =\" l8 b' C( _  m1 ?6 a+ u& K7 B
   
8. 
   : ?  n$ f) w. h4 i; i
9.   1.0086e+006
   $ O7 L\" o1 `\" L; g, e8 w
10. : z7 I: V! y2 V  R  a
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ! c& E' S  O* j/ U\\" W  n% P
4. mvar:2 g) q& |3 x' j- n
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 5 C( F2 y; H7 r* B* j
8. t=clock(),
   
9. & C% ^) r( w$ I; N) J7 y  Q9 }; a+ E
10. oo{\\" s( s4 j4 J, p4 z\\" B1 p% v3 I
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 2 h% |# ^  O1 N1 I6 g, x: Z8 k1 K
14.     arrayfun[HFor("f"),x,y].Sum[]: j0 K- \' L' c0 U* X
    
16. };\\" Q) L7 y4 A6 D$ [, h+ Q! Z
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">! c& E' S  O* j/ U\\" W  n% P</span> mvar:<font class=\\"jammer\\">2 g) q& |3 x' j- n</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">5 C( F2 y; H7 r* B* j</span> t=clock(),<br /> <span style=\\"display:none\\">& C% ^) r( w$ I; N) J7 y  Q9 }; a+ E</span> oo{<font class=\\"jammer\\">\\" s( s4 j4 J, p4 z\\" B1 p% v3 I</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">2 h% |# ^  O1 N1 I6 g, x: Z8 k1 K</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">: j0 K- \' L' c0 U* X</font><br /> };<font class=\\"jammer\\">\\" Q) L7 y4 A6 D$ [, h+ Q! Z</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒

( _# s% i0 d+ z" z* c' v可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
6 u' ]1 Q. u& \( Z3 K
－－－－－－－－

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

我就在你背后

学习中。。。。。。。。。

海水

看看看。。。。。。。。。。。。。

sxjm567

一起交流！楼主给咱们提供机会了