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

forcal

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

1. //用C++代码描述为：* d* C8 Q, l$ D6 V# q) r
   
2. s=0.0;
   2 m8 [# k! @) E* J
3. for(x=0.0;x<=1.0;x=x+0.0011)
   $ \; F# c0 Q3 m. @+ b
4. {' c\" o, D/ Q/ y& f/ G' g2 h
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   , D  o- J. \  t7 R& y, E
6.    {
   3 x: k, A( ^$ A; l$ |
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   7 j+ |& n. u9 f- d
8.    }3 a# k7 r; e) n4 h* r
   
9. }

复制代码

Matlab代码：

1. tic5 [! i* j  j) j2 [4 _
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ) z4 q( y( U) @
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))); z2 W5 O: h7 a4 L: [
   
4. toc
   # X6 S/ {/ g3 X# t1 b
5. ) C; ^& O. T) I) r\" t
   
6. s =
   \" Y# q8 s7 r: M* ~' a' Z% _
7. 
   + g; Y+ x  C1 B3 s  l6 O6 {0 L
8.   1.0086e+0060 C6 v6 R$ z6 N: E  Q
   
9. % x! E5 }, O  L\" c# M\" A+ e2 `) [
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];4 g0 j9 z! f# o+ l
   
4. mvar:# x5 S, \\\" K: J- |7 z/ I' e& J/ p
   
6. t=clock(),! |, ?  L! x+ P* }( `5 I& l4 o; s5 R
   
8. oo{
   
9. 1 X: v' S# n8 Y. U& E5 J* x
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],# K2 C* `( d; J: O9 b5 S/ d
    
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]4 b2 B6 _+ C8 W6 x  o! A8 j
    
14. };
    
15. 5 H2 r# O/ _/ i: I* f
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">4 g0 j9 z! f# o+ l</font><br /> mvar:<font class=\\"jammer\\"># x5 S, \\\" K: J- |7 z/ I' e& J/ p</font><br /> t=clock(),<font class=\\"jammer\\">! |, ?  L! x+ P* }( `5 I& l4 o; s5 R</font><br /> oo{<br /> <span style=\\"display:none\\">1 X: v' S# n8 Y. U& E5 J* x</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\"># K2 C* `( d; J: O9 b5 S/ d</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\\">4 b2 B6 _+ C8 W6 x  o! A8 j</font><br /> };<br /> <span style=\\"display:none\\">5 H2 r# O/ _/ i: I* f</span> [clock()-t]/1000;

结果：
1008606.64947441
' w4 p. R6 m3 H5 M/ `( S0.641

$ A0 V* M& i; i" S7 sForcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。
. _, D/ R7 U4 a9 D% Z1 J" W
& ?# S$ R5 r% |; vMatlab代码：

1. tic
   & Z9 M% w% G' F2 `' U. ^& M
2. s=0;
   ( j; A4 T5 _. h( X5 O
3. x = 0;
   & V7 `# G+ c8 v$ Q1 r
4. y = 1;# M6 W$ `: ?+ R
   
5. while x<1   
   / d. i' v2 P& x9 k) j8 D  D) T
6.     while y<2;        
   # e, B' [. ]# c- Y) a
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ D( q$ E, Q/ `
   
8.         y = y+0.0009;        0 u  i: n) E# [* @  `$ T
   
9.     end- k6 U8 a* a2 t& d1 W  \
   
10.     x = x+0.0011;
    7 g3 G) ]$ R$ }/ @( ^' c
11.     y = 1;\" r1 f% P0 i/ Z' ]
    
12. end3 s# o8 e! t; \0 h$ e8 e5 D
    
13. s
    + s+ Z+ X5 Y' O( x% o* J9 L7 z3 q( F
14. toc
    % C9 c0 I* ~5 d4 O3 w% D( O
15. 
    \" F% l\" n# \9 k' m- _. ]; Z; Z
16. s =. r0 v- a0 b+ }
    
17. 
    6 ^1 R8 a- O' Z/ L  K  Z, R# q
18.   1.0086e+006+ q0 X. z$ c1 e7 T. ^3 c5 R; @
    
19. 
    # x' w. a& d  {. c; d6 S0 g4 b
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:: i& _% E; \' }% ^2 A
   
2. t=sys::clock();1 W\" |+ R/ H7 {( C
   
3. s=0,x=0,
   1 [3 U/ N* O: V$ b! U
4. while{x<=1, //while循环算法；
   8 N8 a. m\" [& W& n5 u/ ^8 s
5.     y=1, # _6 E& r, E  H3 n$ {
   
6.     while{y<=2,
   , r\" J9 K6 D% A- N
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 9 p# F7 L0 O& N  t: y
   
8.         y=y+0.0009 7 h$ X' V6 d& v
   
9.     }, # e; C  Z: G% v2 S+ Q) B
   
10.     x=x+0.0011
    , z8 }) ?& Y% E  U8 C
11. },
    5 C7 A) H$ M* @' ^2 J
12. s;
    5 l/ G! y1 f/ ?, a2 O
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
1 F6 x3 u3 S" N, q0 U+ y6 S" [/ l+ D
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
+ b- |5 h- M1 B$ S9 D9 R4 U" V
2 f% [( n5 i% D" B, B－－－－－－－

/ q1 L7 L7 b4 I0 d1 ?! AForcal中还有一个函数sum专门进行这种计算：

1. mvar:$ P; D5 n$ y\" E5 J
   
2. t=sys::clock();
   7 Y  e; v3 n( b) }( k1 N
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" _, x6 {$ v2 P; A4 d
4. sum["f",0,1,0.0011 : 1,2,0.0009];\" N, X4 V0 o: E4 N4 d8 {0 ~
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
% [+ v4 W8 P0 v+ m1 O0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
8 u. t4 c5 r1 Z. r" f; a& _
( t( ^+ W* ^. @7 b/ F; Jmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));3 P! g$ T+ W3 R
   
2. tic
   ' K6 d* }; t$ {3 b+ p2 O
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, e) P+ S' e* e/ M. a, H
   
4. sum(sum(arrayfun(f,x,y)))
   \" X$ G' c4 i) B3 V% f
5. toc5 [1 X  ?0 Q$ t. u1 _* g! i\" d
   
6. 3 t1 V+ r# x. z/ ]) `
   
7. ans =5 T2 [7 B$ G! t- S6 j2 C5 y6 J) Q' L
   
8. 
   ! A) @( n, o0 H% G
9.   1.0086e+006* m- [6 q& w) v/ e! s8 e
   
10. 
    - T\" E& _8 f/ a
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 3 B5 f: H: f0 f( E' I9 o- Y4 t4 j* |
4. mvar:
   
5. 0 O- ~3 O( K' T+ l% H! [; q
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 4 W( i2 }5 x% D% {
   
8. t=clock(),8 T0 Y( l/ z% @, t5 R
   
10. oo{! J: T  h% q  l! c\\" H4 u
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],3 T  _! j; F% t# w; E+ ?/ F
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. ' Z# b/ D7 {\\" d2 G9 U7 p
16. };
    
17. 7 f! i. U# ^; u& w: J0 e
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">3 B5 f: H: f0 f( E' I9 o- Y4 t4 j* |</span> mvar:<br /> <span style=\\"display:none\\">0 O- ~3 O( K' T+ l% H! [; q</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">4 W( i2 }5 x% D% {</font><br /> t=clock(),<font class=\\"jammer\\">8 T0 Y( l/ z% @, t5 R</font><br /> oo{<font class=\\"jammer\\">! J: T  h% q  l! c\\" H4 u</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">3 T  _! j; F% t# w; E+ ?/ F</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">' Z# b/ D7 {\\" d2 G9 U7 p</span> };<br /> <span style=\\"display:none\\">7 f! i. U# ^; u& w: J0 e</span> [clock()-t]/1000;

结果：
/ J/ x2 U$ r/ p# T% Y- b1008606.64947441
% X, ~; a3 i& r7 N0.735  秒
, R1 u" H7 a+ ?! F; N/ Z3 L( e$ g% Q* e
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

4 R  |5 L( Y$ J& o－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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