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

forcal

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

1. //用C++代码描述为：
   3 o$ |/ q7 [/ U2 r& @\" t
2. s=0.0;
   - O3 V. N5 ]) w; c  m6 E2 O
3. for(x=0.0;x<=1.0;x=x+0.0011)
   5 G3 ]) {0 G# y( f$ V( b% u
4. {$ H; C  [: S7 Q! a6 k
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)8 A/ @) a7 J4 e# b# H& }) R\" K
   
6.    {
   : M! O+ h0 ~! r0 P) C- u& U& f
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ F7 _; {- j\" n6 d, _/ [% x
   
8.    }6 i- Q( `0 s( e) @/ J+ K\" {: z1 E
   
9. }

复制代码

Matlab代码：

1. tic
   ' P+ q/ Z6 A7 r- X4 E
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);6 L, h. V  k5 L; s/ q, R
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))$ K: j  R' P8 _% ]
   
4. toc
   : {+ ^3 H# C3 Y
5. 
   2 B: H# E) K  C# V+ R
6. s =
   6 e0 @\" A, K4 A. l+ l& U4 u& K& u
7. $ X+ ~- u4 V4 q8 `
   
8.   1.0086e+006
   & T# E4 U+ h$ R* i
9. * W# U. g2 C$ k) z' x
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ' F7 c8 ]. S9 s$ Y; r2 x
4. mvar:! G5 P0 r9 x8 M6 ~. q
   
6. t=clock(),* @0 v2 ^9 r4 |7 @9 h
   
8. oo{
   
9. 6 c2 b1 P- L7 P
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: B9 s\\" H- O6 a# D  q
    
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]
    
13. , C5 g$ E; M$ G1 l
14. };
    
15. 3 j0 a5 @* B1 y% ~
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">' F7 c8 ]. S9 s$ Y; r2 x</span> mvar:<font class=\\"jammer\\">! G5 P0 r9 x8 M6 ~. q</font><br /> t=clock(),<font class=\\"jammer\\">* @0 v2 ^9 r4 |7 @9 h</font><br /> oo{<br /> <span style=\\"display:none\\">6 c2 b1 P- L7 P</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: B9 s\\" H- O6 a# D  q</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]<br /> <span style=\\"display:none\\">, C5 g$ E; M$ G1 l</span> };<br /> <span style=\\"display:none\\">3 j0 a5 @* B1 y% ~</span> [clock()-t]/1000;

结果：
1008606.64947441
! O) Y3 ~7 c& p* Q( p# z' v% X0.641
9 l6 n. Q; c1 y/ }
Forcal比Matlab稍慢些。

" t/ ~/ A! n# P－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic: ~  g% X& L. V( R% N
   
2. s=0;
   $ x$ Y: G  x\" K- m1 H
3. x = 0;
   1 z  J3 `/ a# S$ ^3 m! u
4. y = 1;6 x4 Z2 m: i8 X
   
5. while x<1    3 i* t: T  M8 v: O
   
6.     while y<2;        
   ( F) n2 g. I+ k3 S6 \. }0 F, k
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   9 \7 u( G& r. J  F1 W) Z9 X
8.         y = y+0.0009;        ) k3 t, }* y. T8 w3 q
   
9.     end) B4 }: O+ r9 q
   
10.     x = x+0.0011;
    ! z% V( A6 l1 w& W& K
11.     y = 1;
    / M- D8 t. `. ?+ p5 V- H\" m
12. end# I8 c4 A: J+ F9 u9 o- _. c
    
13. s
      l, t& @. G6 C7 |% L, X* ]
14. toc
    $ g! ?/ B* v/ l/ [: W1 a/ d) A
15. - J9 b( {$ d5 J
    
16. s =% D. h8 n$ p3 A) _% r( `- L( t
    
17. 
    ( Q, ^) g3 _$ k& G, z5 f3 u* j/ u; H
18.   1.0086e+006
    \" t8 }; F6 I! m\" f( P& ]
19. 
    & H  s; {- H- t1 N
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:7 a. w0 c8 i0 \, Y
   
2. t=sys::clock();
   / M- k6 V8 t/ ~- e% f
3. s=0,x=0,
   : Y  u! _! N. U0 g+ d1 d2 d
4. while{x<=1, //while循环算法； 3 F+ F' Q1 S- V
   
5.     y=1,
   2 ?8 r\" c: _6 e! l
6.     while{y<=2, - c7 d, b6 S$ y% w\" ?\" X: e1 P. o0 F
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
     O' t: y4 {# Q! y9 @' B
8.         y=y+0.0009 , t% a8 H\" T- _% n4 O4 ?5 s
   
9.     }, 9 u; u. P\" w1 M
   
10.     x=x+0.0011 ) T# k1 N( S4 q  `# D- ?$ p
    
11. },
    + B1 F\" {8 a: v' f: K
12. s;
    7 F: {. n7 j% f+ P( y6 f
13. [sys::clock()-t]/1000;

复制代码

结果：
& C4 U' o5 i5 X% k2 X( S1008606.64947441
2 V0 M& ?4 ~! j3 P0.734   //时间，秒

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

+ x: i4 J( k+ W8 g3 Z! k－－－－－－－
7 h; x- j7 ^' m& W, v+ ^& o
3 ]* Y# ]$ _( C% Y, r/ G0 lForcal中还有一个函数sum专门进行这种计算：

1. mvar:% n$ k( w; B1 \* n8 w
   
2. t=sys::clock();
   0 x) D- j7 W' ]% ^; [
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ; o\" e3 h, e' E- V/ k+ |
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   8 p: B$ B$ \- ?# V: X8 N
5. [sys::clock()-t]/1000;

复制代码

结果：
  p" e9 o! Z2 w4 ^- G9 O: \3 }* D1008606.64947441
6 z2 {( O8 `/ a0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! {6 T8 x9 H$ ]6 c: p, p
   
2. tic; Y( n. n. K  I, Z1 E! f* C
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);' h7 p; Y8 w7 c
   
4. sum(sum(arrayfun(f,x,y)))
   8 q3 k! f/ S1 y! O0 e
5. toc: E) c3 k* D$ _0 x
   
6. 
   * F) G- }) x. M, Y$ R1 d
7. ans =* G; I# ?# @# \* d\" ?( q/ b
   
8. 
   - |6 g1 F  g9 Y: Y
9.   1.0086e+006# @3 w0 P# |( _) A8 n
   
10. 
    6 H9 k: y1 ?0 ~) x) F' z
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. * H! t% h8 \7 q, }0 n! ^$ D
4. mvar:
   
5. $ B8 G' _1 K6 j. o7 g: Y
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); , \* P/ d8 F! A
   
8. t=clock(),& t, ~! o( S/ Y7 M0 |# u
   
10. oo{6 }+ H$ {, r% }* A4 l6 Y
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. & k+ r' H% R% q  i( q% I
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. / D3 c; {3 p; _8 G% `( F% o* C
16. };
    
17. / m5 |6 z2 ~! e& f\\" Y$ O0 r
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">* H! t% h8 \7 q, }0 n! ^$ D</span> mvar:<br /> <span style=\\"display:none\\">$ B8 G' _1 K6 j. o7 g: Y</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">, \* P/ d8 F! A</font><br /> t=clock(),<font class=\\"jammer\\">& t, ~! o( S/ Y7 M0 |# u</font><br /> oo{<font class=\\"jammer\\">6 }+ H$ {, r% }* A4 l6 Y</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">& k+ r' H% R% q  i( q% I</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">/ D3 c; {3 p; _8 G% `( F% o* C</span> };<br /> <span style=\\"display:none\\">/ m5 |6 z2 ~! e& f\\" Y$ O0 r</span> [clock()-t]/1000;

结果：
; p8 Q5 m5 Y+ J8 T7 z6 W1008606.64947441
0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
5 U6 Y+ Q/ |0 p, t$ X% U, p& y
7 e1 m7 N4 @  {% L4 m. O( r% q－－－－－－－－

- i# v9 i2 G; u; `4 T% @从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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