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

forcal

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

1. //用C++代码描述为：% Y! \0 N* d7 d1 j) u  y& f
   
2. s=0.0;
   4 s) q4 T% N  }% n4 c, W
3. for(x=0.0;x<=1.0;x=x+0.0011) 1 |) L8 P2 N( Q
   
4. {
   % h: r0 K1 |1 J9 J
5.    for(y=1.0;y<=2.0;y=y+0.0009)' u) J0 L1 O4 U+ q9 A$ \
   
6.    {$ V' \6 v! A8 J
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));# t' f: }6 P$ m* o. z
   
8.    }
   ) L/ H\" z7 l$ {  t$ w, h
9. }

复制代码

Matlab代码：

1. tic
   ( i7 b' c9 X% _! ?# m
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   7 W2 y& o1 \: L1 S* q+ K
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))): n4 \4 n, |% G
   
4. toc
   ' ^0 j& ]$ x8 m2 X2 a
5. 
   % {) K' p* [, D$ V2 I: I& {. b
6. s =* Q$ G9 V* Y0 B7 i2 \
   
7. 
     f7 Q( H2 G# x5 A! Y
8.   1.0086e+006
   5 H( R0 h7 n1 O  u3 u
9. 
   4 ^0 W- G' Y, T
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];3 w\\" i8 q$ B, x, ~6 j0 q. ]
   
4. mvar:. v% X( A1 U$ R
   
6. t=clock(),$ @8 v; v# _& U1 I
   
8. oo{
   
9. 0 F0 V( C6 h: f1 y+ [3 I  B/ n1 U
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( F5 f& _- j! ~  \# m# d1 N, T0 v( B
    
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]7 J5 Z5 l5 H* H5 `3 s$ q4 v
    
14. };
    
15. 4 U! M/ \, R, Q  z1 W' Q' O$ m
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">3 w\\" i8 q$ B, x, ~6 j0 q. ]</font><br /> mvar:<font class=\\"jammer\\">. v% X( A1 U$ R</font><br /> t=clock(),<font class=\\"jammer\\">$ @8 v; v# _& U1 I</font><br /> oo{<br /> <span style=\\"display:none\\">0 F0 V( C6 h: f1 y+ [3 I  B/ n1 U</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( F5 f& _- j! ~  \# m# d1 N, T0 v( B</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]<font class=\\"jammer\\">7 J5 Z5 l5 H* H5 `3 s$ q4 v</font><br /> };<br /> <span style=\\"display:none\\">4 U! M/ \, R, Q  z1 W' Q' O$ m</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－

1 W$ E5 F9 U# u! C% ~; j6 s再看循环效率。

$ t# c% X# E2 P6 Y' YMatlab代码：

1. tic
   9 w' H$ m3 V# @' {\" {9 Y* A6 I6 J
2. s=0;! I$ k. G% `9 [* U+ n* m4 l
   
3. x = 0;
   0 H' X+ }+ X4 y, f6 c% y
4. y = 1;
   , _2 D# u5 D; W1 a
5. while x<1    2 k& A9 D  S! L0 s
   
6.     while y<2;        
   ! w7 M$ v/ G\" Z0 E; v
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" r* B5 m: V0 p6 i2 {1 a* t4 N4 [6 A
8.         y = y+0.0009;        , f( K6 t% s\" T' K, v, |* s
   
9.     end
   : u4 k3 a( E/ c  Z) A! V
10.     x = x+0.0011;
    ' D/ C/ t8 B' q' E9 }# T
11.     y = 1;
    , r& z# M3 f4 ?+ b/ G; V2 a8 x' f
12. end6 J( `1 i3 }' r: M
    
13. s0 p% N( }) G+ K( }; G
    
14. toc
    : [/ C; {. e5 }- b! |
15. / G  H! ~% Y  W; I\" u
    
16. s =
    & N: J1 C. f* F$ @2 {, x4 d7 q$ P
17. 2 d7 A1 [& _, d
    
18.   1.0086e+006
    ; e+ Q) j# v# @; Y) K7 ?  @' q
19. ) m9 N& H1 A5 C* g6 D  ^6 W5 d\" q
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   ! h+ q$ p) s3 Y7 K7 B' h( w7 o5 r
2. t=sys::clock();
   \" _+ k& C6 |. g/ a8 @* u
3. s=0,x=0,
   ( u7 g& k6 |2 V
4. while{x<=1, //while循环算法；
   * `6 ]* e& F/ T1 V: t. p2 I+ ^+ O
5.     y=1,
     M; H2 v% a/ ]3 S$ V
6.     while{y<=2, & a% I6 u2 t( C& Y% Q2 W
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), * k# n8 I8 O: q# C, |
   
8.         y=y+0.0009 8 I\" |2 l& M1 K3 o\" Y( O# T
   
9.     },
   # p2 _4 a- B- A6 U3 z
10.     x=x+0.0011
    % k( i( P1 R9 y2 L
11. }, & R6 D& B. O% b9 M9 F\" v
    
12. s;- a& x+ b) `+ I! p# w
    
13. [sys::clock()-t]/1000;

复制代码

结果：
5 U! z( v0 T( {: @" |1008606.64947441
0.734   //时间，秒

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

－－－－－－－

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

1. mvar:. \- X: K2 s' N
   
2. t=sys::clock();
   7 b' B7 ?4 U* S( W1 B
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" L+ [2 |- s$ Q0 F4 ~9 Y. K* g4 g  k
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   4 w6 m% n' o9 b/ V
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
+ k" d! s1 E  Q0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) W  n3 b) t( n' n( b
   
2. tic+ ], b8 M6 y6 [8 C; v
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ; p; G1 p3 H& t\" O5 T3 E. [! {0 q
4. sum(sum(arrayfun(f,x,y)))
   0 ^/ u' W* }! Q9 E\" U& Y
5. toc
   2 ?\" R* E: z0 O' T2 h4 G
6. 
   3 V+ l2 p$ @8 b
7. ans =
   ) R0 N2 V9 i- U/ ~2 O
8. % p3 K$ B8 E5 v
   
9.   1.0086e+006
   - |. F% T, v6 g. `& w9 L\" s
10. 
    5 U$ _3 X: X  V0 k( x! q/ L
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];* b4 \7 p: i) F: o
   
4. mvar:
   
5. \\" u  h$ P  j3 `
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 3 c8 ~7 k+ z  k
8. t=clock(),& z) a\\" ^. L: Y6 ^  N) r4 w
   
10. oo{$ K: E$ c) |2 N\\" o' O+ O
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],$ a\\" j: c3 n$ F7 ]1 P
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 3 M; E6 I' m; }, m3 v, x# D
16. };+ c' k3 x% X/ N5 c* p% D
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">* b4 \7 p: i) F: o</font><br /> mvar:<br /> <span style=\\"display:none\\">\\" u  h$ P  j3 `</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">3 c8 ~7 k+ z  k</span> t=clock(),<font class=\\"jammer\\">& z) a\\" ^. L: Y6 ^  N) r4 w</font><br /> oo{<font class=\\"jammer\\">$ K: E$ c) |2 N\\" o' O+ O</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">$ a\\" j: c3 n$ F7 ]1 P</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">3 M; E6 I' m; }, m3 v, x# D</span> };<font class=\\"jammer\\">+ c' k3 x% X/ N5 c* p% D</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
. K6 ?- s# ?/ D/ {. P& O
; i* w( ]  e0 f8 T2 t6 n( B可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

: R5 z8 _! O1 N& J' C+ m9 x－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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