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

forcal

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

1. //用C++代码描述为：
   4 F( O- M, ?# `! v\" q- q0 m
2. s=0.0; ; r/ r) ~6 @( g' `2 X  ^, M! H
   
3. for(x=0.0;x<=1.0;x=x+0.0011) ' ?# s0 ^: ~: D( d# i' t0 D  F; p
   
4. {0 n) u. r7 h  l+ _5 T! B2 Y
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   3 ]% H$ ~! Z  N
6.    {
   8 m) L; X+ V. a3 o
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));5 Z+ C  R4 D3 o/ ]! I
   
8.    }
   , s. x$ q$ J& c2 o4 w0 i
9. }

复制代码

Matlab代码：

1. tic
   # j% ~& ?$ K/ x1 |
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);* K* E: I4 t# R* j8 z% s: \
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))9 V7 r% G/ o2 j2 U- I3 ^
   
4. toc
   $ Y1 ]$ w- Z. B* [- L) s; N
5. 0 S* |' i9 H2 S8 |6 D
   
6. s =
   ( |; n4 l: S  A
7. 
   2 o1 ^) c, O  l& ^! n% h
8.   1.0086e+006
   : z, W. c. Z* \$ O0 b4 J
9. $ P7 ^, S5 \% Y1 c
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];( D: ]- P\\" _+ r( R; `; {
   
4. mvar:
   
5. 0 k! a  T0 K7 I4 P1 ^
6. t=clock(),
   
7. , T8 D4 J/ ?+ N. b7 f+ t) C
8. oo{/ O' W+ p+ E/ N1 u* d' D& h7 |
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 6 [7 p; p; c, |) H
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]% X, L6 Q; n; j! D0 J0 A8 D
    
14. };2 T; ~6 s/ j! @5 e# j
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">( D: ]- P\\" _+ r( R; `; {</font><br /> mvar:<br /> <span style=\\"display:none\\">0 k! a  T0 K7 I4 P1 ^</span> t=clock(),<br /> <span style=\\"display:none\\">, T8 D4 J/ ?+ N. b7 f+ t) C</span> oo{<font class=\\"jammer\\">/ O' W+ p+ E/ N1 u* d' D& h7 |</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">6 [7 p; p; c, |) H</span>     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\\">% X, L6 Q; n; j! D0 J0 A8 D</font><br /> };<font class=\\"jammer\\">2 T; ~6 s/ j! @5 e# j</font><br /> [clock()-t]/1000;

结果：
6 S/ E0 a- u; a7 w2 r1008606.64947441
0.641
  W* C7 q  q2 q' v& @* g, E
Forcal比Matlab稍慢些。
4 d/ j8 L% @4 E: W7 c2 m1 w
－－－－－－－－－－
: c/ ?8 E5 P, P# R  L0 R: L6 A' e
1 z: ~; l  ]; j. \再看循环效率。

* {! R. O& h9 Q( W" g' S- EMatlab代码：

1. tic
   2 k2 M+ R2 m6 W' C# l' N
2. s=0;, H7 G* c0 \  Z9 b  }
   
3. x = 0;8 j% W5 u  V& L) @5 q
   
4. y = 1;6 c, d2 ~3 F! v& ?4 j- M
   
5. while x<1    4 H) o* I* m( M# U1 W2 y
   
6.     while y<2;        
   # J: b/ A: d1 P$ q; W8 G# j% O
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! ~: }: K- w$ J; f+ E& v; U
8.         y = y+0.0009;        
   9 J% _, I, J( h8 Q+ _
9.     end/ B, B4 r# C, V0 J4 G& x8 d; Z2 u
   
10.     x = x+0.0011;6 E6 }\" |; t\" m7 T% q; N
    
11.     y = 1;/ j% ^1 ?- }4 k9 g
    
12. end3 t% J; K3 g+ u, H1 @$ ~\" B8 c
    
13. s
    . C: ?# c% d: Q+ U- G
14. toc
    ) f/ v% A! D. y3 S( W
15. 
    8 r2 [& M  G% u6 h9 R$ [
16. s =
    \" U) y: a% u( _1 ^6 A( z
17. 
    , K7 N5 e4 `+ s
18.   1.0086e+006( o+ v5 y2 D- \: U6 `  X
    
19. : R# h$ ]; C4 U6 {* [
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:3 k\" D1 d/ W5 h& G, ^/ b. [& g
   
2. t=sys::clock();
   5 Z2 ?+ S  [! s) Q) g\" e
3. s=0,x=0, # P  v; R. D$ ~2 {5 Q# w
   
4. while{x<=1, //while循环算法；
   $ b, ~\" @, i* Y
5.     y=1, % c+ ^( O4 z- u0 F* }
   
6.     while{y<=2, + I$ K2 P7 u. d8 Q
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   $ N0 F- |# V: n' i' O
8.         y=y+0.0009
   6 s% u8 ]& ]2 R
9.     }, ! D: ~; m& u' z/ [: s' i
   
10.     x=x+0.0011
      s\" @) S$ }. q
11. },
    6 O' l- F* N/ e! A
12. s;3 k$ j1 b/ n. u' n) m' y3 f5 U  t
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
$ t9 B; z8 O* |7 J: N5 s! B' G: i' W0.734   //时间，秒
2 X6 `) U0 Q! q  M/ g; c) Z9 h% i
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

5 R/ \+ G2 I/ [& ]' n- H; g5 T－－－－－－－
  n/ }/ T4 d  g1 w0 ^, E# j
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:7 t% V9 ?, C6 y( i- h$ M
   
2. t=sys::clock();
   # Z1 h5 I2 I( [1 u. j9 t5 m\" N
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . m, t+ w5 l9 h/ q
4. sum["f",0,1,0.0011 : 1,2,0.0009];9 @% b\" |\" c/ s# p1 L  R: P2 U
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ u2 v1 p# d4 ~: Y7 k
   
2. tic! j3 a' Z- Z  u. f4 |( ]9 O
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   - d0 q! _$ W6 n0 x0 ~
4. sum(sum(arrayfun(f,x,y)))
   % E3 j- V, l8 L: J; s. `1 ~
5. toc- C2 W- f, e4 G
   
6. & @. B- o9 c1 m/ F+ \5 a2 Y
   
7. ans =
   + ]  p! [. r  }\" V* C
8. 
   # [: Q9 a' t- E' a+ ~
9.   1.0086e+006: j# `2 b6 ^0 P+ q
   
10. / U% [. c4 r3 q\" l  F& g
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];3 g5 W( G  ]; g8 ^5 V
   
4. mvar:
   
5. , h# F- g' n9 P( f
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); / Z  e! c- i7 q) A0 k\\" U& v
   
8. t=clock(),
   
9. % V: w) h; Z( b1 r% z) \
10. oo{
    
11. 5 U. V6 ^( B1 C+ v. m
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 6 g( P8 D, ^, x, X: H4 N8 p/ |
14.     arrayfun[HFor("f"),x,y].Sum[]+ s4 _- C: u9 c4 s1 W\\" w3 q1 U
    
16. };$ H% n2 t0 b* S0 A2 m
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">3 g5 W( G  ]; g8 ^5 V</font><br /> mvar:<br /> <span style=\\"display:none\\">, h# F- g' n9 P( f</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">/ Z  e! c- i7 q) A0 k\\" U& v</font><br /> t=clock(),<br /> <span style=\\"display:none\\">% V: w) h; Z( b1 r% z) \</span> oo{<br /> <span style=\\"display:none\\">5 U. V6 ^( B1 C+ v. m</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">6 g( P8 D, ^, x, X: H4 N8 p/ |</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">+ s4 _- C: u9 c4 s1 W\\" w3 q1 U</font><br /> };<font class=\\"jammer\\">$ H% n2 t0 b* S0 A2 m</font><br /> [clock()-t]/1000;

结果：
3 \. c5 r! Q- @% o+ k$ B; f- U1008606.64947441
0.735  秒
* W6 \: q/ {0 r% ?% ~
1 h- f/ {: p: j3 Z! M. J+ f可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
2 H/ N! U  P. n* m
－－－－－－－－
% U" \0 o3 b) N7 g
" A6 R7 m9 ~& l* D% j8 f从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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