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

forcal

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

1. //用C++代码描述为：
   * d( v5 S7 d\" v6 U
2. s=0.0; 2 B* V+ ?\" t2 r* h( @: i5 y: q
   
3. for(x=0.0;x<=1.0;x=x+0.0011) 6 h; b6 T* m: j( a/ k3 H
   
4. {\" U5 @$ p: P- s4 [: Q8 J
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)/ D; y7 {* K0 r1 ]( w
   
6.    {
   + Y, U2 z2 J1 P9 y* j8 J+ |
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 m+ |5 ^5 L* n+ |, u* D3 Y
8.    }
   / w# D3 ~  T% F6 M6 O* G
9. }

复制代码

Matlab代码：

1. tic/ x  h, c0 O# ]( {0 v9 d9 B
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * Y% q, E5 e, P; |  O* ]
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))/ e1 x) b4 [; s  b0 V/ x9 |! Y
   
4. toc3 L, \! e0 ?5 s1 F7 P
   
5. ' }# {- z/ n) ~\" i3 M
   
6. s =* A+ i! u8 a. e' C4 [: E! W( M
   
7. 
   9 }) ]8 G\" t# D/ H5 Z8 t0 c/ o
8.   1.0086e+006; T9 n, }) v$ m  n9 \' x
   
9. / t' A* I6 ]5 S, A: c
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];5 q) Y: N2 K2 L' U
   
4. mvar:/ L! ~3 |\\" B: P# g6 `: v# N$ x( K
   
6. t=clock(),
   
7. $ X! Y/ N$ U) a# p' `7 j\\" z
8. oo{7 i# r) y' j0 d+ E. G9 l
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],/ n0 t: @( q' f& 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. 3 e( X8 F$ r; W& |9 B' o9 y# ~
14. };5 E$ b, K; o; @1 O, ]7 J, e4 R3 k
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">5 q) Y: N2 K2 L' U</font><br /> mvar:<font class=\\"jammer\\">/ L! ~3 |\\" B: P# g6 `: v# N$ x( K</font><br /> t=clock(),<br /> <span style=\\"display:none\\">$ X! Y/ N$ U) a# p' `7 j\\" z</span> oo{<font class=\\"jammer\\">7 i# r) y' j0 d+ E. G9 l</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">/ n0 t: @( q' f& 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\\">3 e( X8 F$ r; W& |9 B' o9 y# ~</span> };<font class=\\"jammer\\">5 E$ b, K; o; @1 O, ]7 J, e4 R3 k</font><br /> [clock()-t]/1000;

结果：
3 T; y; p8 f+ r- I1008606.64947441
) T# H8 K2 x) u3 p) n7 S0.641

Forcal比Matlab稍慢些。

2 K3 v0 |$ @' C－－－－－－－－－－

& \0 b; m& e, f3 S8 f2 N再看循环效率。
  R! u: {& a) R3 @$ C) {
2 b& Y. `/ k; o# p9 [2 Q9 u' \5 hMatlab代码：

1. tic) p7 n( p6 H7 A
   
2. s=0;
   8 t: j6 d/ \% l5 P1 d, ?
3. x = 0;! ^& w7 o7 P& j) C4 \+ w, |
   
4. y = 1;
   9 \# w\" Z$ _- |3 r2 y
5. while x<1   
   3 e* Z5 |7 l1 k1 [) |
6.     while y<2;        
   0 y% x/ t* y2 i! q
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));. l; S& E  |; I9 I! u, z
   
8.         y = y+0.0009;        ! C8 K; ^5 W* _/ l
   
9.     end
   , M6 |9 y5 L: }
10.     x = x+0.0011;% L: i# L- M/ w! j5 u! H
    
11.     y = 1;
    ( C6 c9 R+ Y- n
12. end: \( V8 k. J& e4 A9 t) s% G
    
13. s
    5 F% s\" ^/ u( W* b5 Z
14. toc+ C6 c. ]+ _, m0 G0 G* `
    
15. ; L( b( M4 V2 k6 @
    
16. s =9 g1 W7 {) m: v7 Y+ U
    
17. 6 h  y9 N2 ^1 n/ L; V. B
    
18.   1.0086e+006
    , L, B- v. d( p+ ^( F
19. 
    4 @! x; z; ]  p$ R: e& p: ]. g9 t
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   , i+ k1 `, N! k
2. t=sys::clock();
   - O0 A  k0 ?  A' z- p# G
3. s=0,x=0,
   * E0 y& |; m  ]\" b8 P, V
4. while{x<=1, //while循环算法； * ?' }, s9 C+ n5 W$ M# D8 L
   
5.     y=1,
   / j9 A4 Q6 D+ g7 |5 g1 M3 g& {
6.     while{y<=2,   I# e: m/ q9 O. {\" P* G
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 0 ~+ K\" l, f9 e7 G- u
   
8.         y=y+0.0009
   ! ?# C+ D! J9 Z2 ^. y* @0 g( P
9.     }, . b3 N6 x1 n# [
   
10.     x=x+0.0011
      g' I$ n' i! v3 D$ c' s; r5 g
11. }, 1 y! [4 I5 F5 {\" r
    
12. s;
    ) b. @6 \! P  B\" e  o, ]5 N
13. [sys::clock()-t]/1000;

复制代码

结果：
7 `( G# F: i( `  M, j1 W) K2 y1008606.64947441
0.734   //时间，秒
% H# w, I) N  b2 X$ \8 J( Q, m
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

* ]7 j' |2 t& \& x9 z, ~+ }－－－－－－－
) X/ Q; X+ g, c' o3 l5 x( W) y
" |: W: R: Z8 e$ \Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   & y2 q/ x\" B  p
2. t=sys::clock();- a6 w8 L: J$ D$ q9 a. G; W
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 3 |* c% I. p7 Q$ ^% r7 r4 A\" m
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];/ J* t6 |- ]2 B1 w
   
5. [sys::clock()-t]/1000;

复制代码

结果：
2 u6 L* d2 c( ?) L7 T- U1008606.64947441
0.719   //时间，秒

forcal

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

  M6 W  e2 b7 U$ jmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  f7 j- y+ E! e
   
2. tic
   $ F. ]: K$ O' M1 X2 x\" }8 K' V
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   # s+ o5 ~5 U/ Z6 u
4. sum(sum(arrayfun(f,x,y)))\" `' t# c) ^1 Y
   
5. toc  |' B3 A# N5 S3 l+ ?# v
   
6. 
   9 `4 x) f/ a! _0 L+ ]! z; d
7. ans =
   : Q\" s' F7 w% i2 A
8. 
     G: v4 l; a$ A+ ?; a
9.   1.0086e+006
   5 }% n( I3 d4 z
10. 
    1 V\" |2 B; I6 P2 p1 a* Q& B( X
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];8 C# ~1 h) r- O( l9 }& C4 a: ~
   
4. mvar:- l$ j: U5 s1 n
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. , C' w9 }7 b4 \( ~: W\\" L
8. t=clock(),7 C  y2 ^, ]% y. e: }
   
10. oo{, f. i- T\\" r( H$ h  X
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. * L+ K; t& y  E4 A
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 3 v7 {4 b. y$ [9 A1 p
16. };  y, s1 Z9 v# U0 O
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">8 C# ~1 h) r- O( l9 }& C4 a: ~</font><br /> mvar:<font class=\\"jammer\\">- l$ j: U5 s1 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\\">, C' w9 }7 b4 \( ~: W\\" L</span> t=clock(),<font class=\\"jammer\\">7 C  y2 ^, ]% y. e: }</font><br /> oo{<font class=\\"jammer\\">, f. i- T\\" r( H$ h  X</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">* L+ K; t& y  E4 A</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">3 v7 {4 b. y$ [9 A1 p</span> };<font class=\\"jammer\\">  y, s1 Z9 v# U0 O</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
7 E& D( ]+ t9 J1 P/ l5 u
! w4 g: Q% z9 B8 h% }可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
, |1 a' H) ~; S* s
－－－－－－－－

- ]6 r, o9 H- t4 ]. d& X$ Y) u7 P从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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