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

forcal

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

1. //用C++代码描述为：
   \" H# N- `3 }/ p  }) m8 R( u2 q
2. s=0.0; ; H# ]2 ^- K) C$ v+ D. L  h3 n
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   . k7 \/ {2 p( U\" W( w) b- {1 O- Z- y
4. {
   % ^5 @. i) L$ D  c+ @$ J& V3 b7 h! R8 C
5.    for(y=1.0;y<=2.0;y=y+0.0009)) o0 j- Q0 y( p4 r; ?. @\" _- {
   
6.    {3 r; d2 e3 x9 m) m$ y' u
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) Y, ^) I# B3 k
   
8.    }, w# s; g: k+ k3 G
   
9. }

复制代码

Matlab代码：

1. tic
   5 n: ?: W  B3 h' _1 I7 R8 D
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);9 W4 _. k# P/ T
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))/ z0 t1 f! Q/ B( o( {, Y; n, c
   
4. toc\" W5 Y7 l# t\" ~+ q/ I
   
5. \" Y6 G) t+ m\" x; F( n\" R( W/ Y
   
6. s =4 ]  @* S- p' U' R$ ~
   
7. 
   , l2 E$ d( h$ l* b  v7 i; N/ |
8.   1.0086e+006. Q6 z2 w: X. n5 R+ j! S- l+ Y
   
9. 
   + I1 x4 b/ r& j( x
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ( c0 O6 {3 C* j& Q- c4 |9 Q
4. mvar:
   
5. & {- J% ^* ~, y( j! z
6. t=clock(),4 [( V: E  M* w2 b( I
   
8. oo{6 c8 ^- `% T  q
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: M0 a; N  B+ Z# J\\" h6 P
    
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]( M% U' M5 _! t2 R0 z
    
14. };& F( k: e: L' n5 V; |- O  [
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">( c0 O6 {3 C* j& Q- c4 |9 Q</span> mvar:<br /> <span style=\\"display:none\\">& {- J% ^* ~, y( j! z</span> t=clock(),<font class=\\"jammer\\">4 [( V: E  M* w2 b( I</font><br /> oo{<font class=\\"jammer\\">6 c8 ^- `% T  q</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: M0 a; N  B+ Z# J\\" h6 P</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\\">( M% U' M5 _! t2 R0 z</font><br /> };<font class=\\"jammer\\">& F( k: e: L' n5 V; |- O  [</font><br /> [clock()-t]/1000;

结果：
4 e0 q  ], s% t& o+ k8 j, B1008606.64947441
- l% \/ c6 T2 E5 |! `6 [; S0 ?" m0.641

Forcal比Matlab稍慢些。
# P) e% A/ e; h2 g, }
－－－－－－－－－－
" Q+ _5 ]8 y: b
再看循环效率。
& z) p: a) Q8 v+ v! X: l4 g8 s
Matlab代码：

1. tic
   + y+ Y6 T% p$ c2 e& r
2. s=0;% Z3 f8 J$ Y2 Y$ n
   
3. x = 0;
   , g\" r+ y1 d6 ?( B0 N/ A- H
4. y = 1;
   / |; F4 R. ?3 e
5. while x<1   
   4 H\" r; S0 W- M7 b- }
6.     while y<2;        
   $ G: [# ]5 r6 }* s
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & L, W# ~, |$ P
8.         y = y+0.0009;        ' T! L$ R* A7 ^) R3 L& |! i
   
9.     end# n$ ~0 B( J/ G6 A
   
10.     x = x+0.0011;
    . U\" W4 c8 Y6 [
11.     y = 1;, A) ^- Y  i& ?4 {! F4 _0 B% o
    
12. end
    0 f8 _' F\" b- ~4 d. Q, C
13. s1 f5 m2 t( |! M
    
14. toc
    4 h+ W: g/ [: H$ _8 o0 c$ p$ s3 y/ @
15. 
    ; ]7 ?' \, b. l\" }7 b
16. s =  k2 Z7 v# N- W; T
    
17. 
    7 D2 ]8 l4 h0 r\" k# x
18.   1.0086e+006
    , Q* C% K! b$ D\" v
19. 
    0 c9 ~6 o, _9 V& Y
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   1 s! X  q& `0 I1 C\" ~; Y
2. t=sys::clock();
   4 d$ b- v5 g! t
3. s=0,x=0, ! ^' R: ?  q6 [9 S
   
4. while{x<=1, //while循环算法； 9 E- R* m( ?4 c  y) X
   
5.     y=1, : x2 u# y. A, _% D  Q4 s5 n. f
   
6.     while{y<=2, 8 P( D3 u! v) Q( C\" y: y
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),   M  T0 N4 ?\" X4 c0 ~7 T% R% |8 o
   
8.         y=y+0.0009
   * n; u1 P4 @: @% |
9.     },
   ' R6 p; q& w$ S/ k
10.     x=x+0.0011 ; H0 ~0 i8 ~* }\" E
    
11. }, ) y, U- @$ x3 y* b/ \
    
12. s;3 U5 j  L$ q\" w7 ?\" ~
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
" Z3 H8 P$ e' v$ n1 f: @, m6 L2 S0.734   //时间，秒
  j+ ^9 \$ x/ C* z( Y5 R3 K
$ o& f! d6 ?# ]- J  ]" c我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
% @- p) Y/ p% w/ d! w/ G
5 B8 g+ w, u' S; H5 xForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   / c) k2 q* k0 |  M
2. t=sys::clock();
   4 k; ^  N5 N: v
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); * h5 a. }# a# z- u+ q\" _
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];. e/ S) R, G\" E
   
5. [sys::clock()-t]/1000;

复制代码

结果：
" Y& @' d6 }# x2 Y$ t. L. u8 S! ~4 C: H1008606.64947441
: V5 k+ t% e9 d: g9 T2 M) [0.719   //时间，秒

forcal

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

1 U, Z1 v  J- rmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));3 H5 ]3 V: l/ F# ?1 F' f
   
2. tic5 ^& H) l5 ?9 P- Q2 ^0 Q- Q
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   6 K* x; n$ m\" i$ b& j( t
4. sum(sum(arrayfun(f,x,y))), z9 [8 T# ?# C1 V' Y! }
   
5. toc
   : q2 O# n\" ?1 Z$ l( H% l; J  K
6. 
   0 p1 R' d# p2 ^; ?0 K
7. ans =\" R. d. E( ^) j' l% ~5 U8 u
   
8. 
   4 m5 b2 x$ P7 v9 J$ Y$ T1 A( Q4 P9 B8 Y
9.   1.0086e+006
   1 ~5 u2 x\" P/ a) @5 V
10. 
    ' q6 A3 I2 H' I2 g
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];\\" U: U8 c6 K7 m! `* V
   
4. mvar:4 `# }* e/ ]% S) F7 F4 }9 W2 B1 c
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 6 A& V+ a8 y2 e& t
8. t=clock(),
   
9. # i7 \\\" w2 c1 x+ A& U( Q3 O
10. oo{7 _; O4 g# z+ a, o\\" @  p
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. $ X2 m/ v6 O( z( V$ p% G\\" ?+ q
14.     arrayfun[HFor("f"),x,y].Sum[]6 r$ w' [8 Q6 a7 S$ M\\" T
    
16. };\\" W8 Z' x7 Y) t
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">\\" U: U8 c6 K7 m! `* V</font><br /> mvar:<font class=\\"jammer\\">4 `# }* e/ ]% S) F7 F4 }9 W2 B1 c</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\\">6 A& V+ a8 y2 e& t</span> t=clock(),<br /> <span style=\\"display:none\\"># i7 \\\" w2 c1 x+ A& U( Q3 O</span> oo{<font class=\\"jammer\\">7 _; O4 g# z+ a, o\\" @  p</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">$ X2 m/ v6 O( z( V$ p% G\\" ?+ q</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">6 r$ w' [8 Q6 a7 S$ M\\" T</font><br /> };<font class=\\"jammer\\">\\" W8 Z' x7 Y) t</font><br /> [clock()-t]/1000;

结果：
2 x) _1 H- H7 u( g+ \8 Z! U# U. P! N1008606.64947441
0.735  秒
3 ^, B$ q2 I  {0 k" n" H) D! U
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
/ C+ `$ b2 \: U1 q
  ]; D) J7 _: {8 a* B从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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