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

forcal

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

1. //用C++代码描述为：
   & U& L* }9 N\" L
2. s=0.0;
   1 Y, {1 k\" L\" P/ ^+ X
3. for(x=0.0;x<=1.0;x=x+0.0011)
   8 }+ h7 `( A2 t
4. {: {% j1 `, d) ~3 u0 i
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   7 U6 j/ _8 g2 [6 F
6.    {
   8 D1 I4 m\" V3 X4 Y! e& Q2 q
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));2 _* g! t0 ?8 Q. @' \
   
8.    }
   % k1 u) @/ K3 m4 F) }$ Q! ?
9. }

复制代码

Matlab代码：

1. tic7 W% D$ g6 z; `  J- T2 N) k' j
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);( r* E5 P6 Q* b% o4 o3 ^& R
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))8 x* f9 i% ^& e5 w$ y
   
4. toc
     l# U1 h\" Q' c, l2 y6 u
5. 
   ) z3 v( W( ?$ w+ @! W; m
6. s =6 D8 `! ]1 ]0 n' U7 y
   
7. * F3 k+ Y$ D* W% [9 a; p
   
8.   1.0086e+006
   \" ~\" l: B6 {) U' O0 ^4 L
9. * g5 ^# F9 G$ t- F1 f8 A
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. + z/ @$ O' p( E8 E, Y% o2 N
4. mvar:4 F! _! v% W! ]; M6 A
   
6. t=clock(),& Y* e: C. |/ k- U7 ?
   
8. oo{7 r. u* q8 e\\" v: X0 ?: c. |9 a1 e
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 4 o0 }' g  R7 k0 I$ e; b4 O- ]
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]9 o, Y$ L! S9 }! k; Z
    
14. };- b7 \! q5 w' \9 {& s; {
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">+ z/ @$ O' p( E8 E, Y% o2 N</span> mvar:<font class=\\"jammer\\">4 F! _! v% W! ]; M6 A</font><br /> t=clock(),<font class=\\"jammer\\">& Y* e: C. |/ k- U7 ?</font><br /> oo{<font class=\\"jammer\\">7 r. u* q8 e\\" v: X0 ?: c. |9 a1 e</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">4 o0 }' g  R7 k0 I$ e; b4 O- ]</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\\">9 o, Y$ L! S9 }! k; Z</font><br /> };<font class=\\"jammer\\">- b7 \! q5 w' \9 {& s; {</font><br /> [clock()-t]/1000;

结果：
  t: c" E- J6 }4 t2 w1008606.64947441
" ~# h$ V& P( l+ ?5 F1 \, q0.641

8 E& Y! H! L# r6 [) o2 }% m- rForcal比Matlab稍慢些。
- ~; \1 `! H; y# X3 b& o
  z) E+ c9 S3 [9 r5 p3 d* k8 g& G－－－－－－－－－－

再看循环效率。
  z( x' Y$ N' G/ d) L+ n
" }  G& b3 J, P9 V8 {Matlab代码：

1. tic# _! F2 N5 V8 W8 g% D
   
2. s=0;; b6 m0 Z9 b* e5 r8 I\" }4 c7 K
   
3. x = 0;
   9 o( ]4 e, K/ ^6 h0 k4 V
4. y = 1;$ ~8 o1 Q3 e8 Z% ]
   
5. while x<1    1 j: O4 E1 H\" ^( A2 k
   
6.     while y<2;        
   ( W, y, ?, B2 r& [- K& B3 y
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / N/ a7 s1 S$ E4 P3 d
8.         y = y+0.0009;          o' N- ?; Y- k
   
9.     end
   , k8 j% C- H& i( E$ f  c5 F
10.     x = x+0.0011;& ]: K\" u& u, ?$ b  G
    
11.     y = 1;
    7 U6 @7 y% C6 @1 v+ Z) j8 \6 {: g
12. end% d/ n# Z0 W3 L
    
13. s
    4 g0 J- t2 y- U- h  Z3 r! h
14. toc
    2 X\" N# {2 P* J/ {) I
15. 
    : C2 f; l1 D0 O2 X9 ^! J
16. s =
    ( F( N6 }' d8 M* w. r' F' J+ c
17. 
    \" ~2 b8 e+ ^( g1 _7 V# q
18.   1.0086e+0063 V2 C: s3 S  _
    
19. 5 i+ |% \1 G1 P* n) K6 k! P
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   1 {6 ]  `& D/ \& a/ V1 W
2. t=sys::clock();
   * C) X* Y( P% F9 ?$ G: I
3. s=0,x=0, 7 q6 k3 l& M\" s
   
4. while{x<=1, //while循环算法；
   ( ^& n! G, u1 v! l& w
5.     y=1, - M! @1 k$ ~! [0 @
   
6.     while{y<=2,
   9 Y, s5 d2 i3 P* ~1 \; H  S
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   * k; I: E/ P6 E  m7 R2 h
8.         y=y+0.0009
   & z, B5 L* t5 m9 m
9.     }, 0 P( n, a\" a) ~- U
   
10.     x=x+0.0011 0 t9 S; R+ W3 |# A& `9 D9 S
    
11. },
    * {2 w: l% w5 K
12. s;  B: ~2 f5 i; W$ o1 ~
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

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

0 L4 F$ L, e/ \- r1 W/ C－－－－－－－
( a* ]: F. n7 f( F& j$ s
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:  s4 _8 [3 J4 F2 y
   
2. t=sys::clock();
   . g: K: @& |' L- j# R4 Z
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); . F9 }# [( C7 n  [
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];; c- W# d- Y\" h+ J. {
   
5. [sys::clock()-t]/1000;

复制代码

结果：
, |3 r/ Q* p  e+ X0 m* _, f! i8 c1008606.64947441
# {7 E$ e% g5 n7 I, J$ {0.719   //时间，秒

forcal

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

5 B% \! X/ h; ?' ~# o2 b) y  Smatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( q( E- j8 Q; W+ e2 |4 z1 ~
2. tic! R9 \$ N\" C: h0 g$ [) W) j- u
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);! \\" `  o! f; K5 S' e, c: Z- T. K
   
4. sum(sum(arrayfun(f,x,y))); S+ I5 X9 x& f
   
5. toc
   3 e, k, B3 A3 s8 ]
6. 
     b7 i0 _: q( F\" @+ ]
7. ans =
   2 {& G+ b' ^* ]( G
8. 
   \" i0 J3 C3 j5 P+ q
9.   1.0086e+006
   ; J3 x7 d; |# K% r
10. : ]/ ]# P) i3 D$ a+ D' X/ ?# t: y. D
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];& b2 c/ C6 `8 R; _
   
4. mvar:
   
5. , T2 N. z4 f9 E
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. : I: E4 j' V$ h7 {6 d& H
8. t=clock(),
   
9. 4 e\\" q: F4 C  N* t, S8 S
10. oo{# x* r) ]3 D7 j7 O9 u) J
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: C, [\\" P# [' N! ]' ~+ |8 A, s
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. $ B- @\\" N* B( g
16. };
    
17. 1 d' B1 _1 w9 p
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">& b2 c/ C6 `8 R; _</font><br /> mvar:<br /> <span style=\\"display:none\\">, T2 N. z4 f9 E</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\\">: I: E4 j' V$ h7 {6 d& H</span> t=clock(),<br /> <span style=\\"display:none\\">4 e\\" q: F4 C  N* t, S8 S</span> oo{<font class=\\"jammer\\"># x* r) ]3 D7 j7 O9 u) J</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: C, [\\" P# [' N! ]' ~+ |8 A, s</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">$ B- @\\" N* B( g</span> };<br /> <span style=\\"display:none\\">1 d' B1 _1 w9 p</span> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
% F  a' J( _) m  v4 B
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
, i( {6 @" u0 p2 q6 \: e/ R/ J+ U
－－－－－－－－
6 Z. x1 v# d" o  |3 }; g& b
* A% {" S8 O! A$ o  i从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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