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

forcal

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

1. //用C++代码描述为：
   ; |/ F4 Y1 o) ?7 B/ M0 B( J
2. s=0.0; + r1 x3 h6 c, I9 {4 T
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   5 ^% Z; L: t% w3 w
4. {
   , L7 Q\" V' f8 e5 ?7 W# b1 S7 @
5.    for(y=1.0;y<=2.0;y=y+0.0009)) _1 B/ p- H+ W8 b# ?% |
   
6.    {! l9 S1 y# [7 |+ w
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 u  q6 @! c. N6 Y\" H
8.    }! P9 D' \- O) U+ `  A$ C0 ^9 e
   
9. }

复制代码

Matlab代码：

1. tic7 N8 }2 v8 H7 L2 i+ S
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   % }9 h. G# c0 ~0 Y; {+ N
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   % g% {/ f2 f4 j: ?  T- i' K
4. toc6 X3 D% G0 Y) u
   
5. % f( ^% N+ r) m  k! m  T
   
6. s =
   0 }) r\" }1 b7 R/ k
7. , }% w2 P) n# G; @
   
8.   1.0086e+006
   % G% F1 ~2 ~$ Y- f
9. 
   1 B1 p1 H$ W+ p0 F: I
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];3 m; u* m: V. t# a
   
4. mvar:
   
5. \\" a# b  I4 E$ ]4 Q
6. t=clock(),8 M/ l% E& |  c( [\\" I6 R
   
8. oo{: i' c* h! x! T8 z# s6 T0 b
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],0 ~& D1 K  h) A
    
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. : H/ M) A4 K8 W4 q) K2 c) l3 {
14. };
    
15. ) K) n6 Z  u0 i9 A8 [8 X3 P( e0 m& ^* Z; `
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">3 m; u* m: V. t# a</font><br /> mvar:<br /> <span style=\\"display:none\\">\\" a# b  I4 E$ ]4 Q</span> t=clock(),<font class=\\"jammer\\">8 M/ l% E& |  c( [\\" I6 R</font><br /> oo{<font class=\\"jammer\\">: i' c* h! x! T8 z# s6 T0 b</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">0 ~& D1 K  h) A</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\\">: H/ M) A4 K8 W4 q) K2 c) l3 {</span> };<br /> <span style=\\"display:none\\">) K) n6 Z  u0 i9 A8 [8 X3 P( e0 m& ^* Z; `</span> [clock()-t]/1000;

结果：
5 R3 n+ u1 D: e5 x* l" @1008606.64947441
' G) G+ Y5 X' {$ p& q2 @: L+ i6 ^0.641
" K9 [0 i9 g) [' w1 F4 [
/ X& E1 n2 R- qForcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic1 M) D! t0 }0 X$ _
   
2. s=0;
   ! S2 O# _. O# O2 \1 `
3. x = 0;! _. Q/ t; P! s2 b2 J\" |
   
4. y = 1;
   $ W% v: V8 `. \* E2 @
5. while x<1   
   & J7 P# s. Y8 S% e8 q
6.     while y<2;        $ _' B& J! B) Z- u1 I9 Z5 N
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ R, ?% I9 i  T\" B
   
8.         y = y+0.0009;        
   # Z- j1 K5 O0 O6 B/ f  W( a8 [- U+ O
9.     end
   . B1 N% D/ }: n1 q/ f% S
10.     x = x+0.0011;
    - _0 U. {! W( q\" w
11.     y = 1;
    : I) G6 C9 g( n1 e4 [
12. end; z1 R4 }( U0 t9 r9 F
    
13. s  B7 T) f& k9 b8 b. g
    
14. toc2 K& B/ X  s; ^4 ]( @
    
15. 
    7 p) p: ~- R# {3 h. z5 b- \' H
16. s =* K6 ]% x+ K/ f
    
17. 
    6 F: _6 k$ v0 c  p
18.   1.0086e+006
    3 F0 F2 |7 {\" h8 A: |4 J
19. 
    ! N  P$ |, P6 {# f9 f
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:3 U# r  m2 ^. K* t: \* }# X
   
2. t=sys::clock();
   7 M$ B8 @, u$ m0 e. H7 x
3. s=0,x=0,
   ) E0 u, S; k( {# S$ J  o
4. while{x<=1, //while循环算法；
   4 B) }; g' G6 |5 E
5.     y=1,
   + _7 Q1 D0 s& a; a
6.     while{y<=2,
   ! C. Y7 M0 H\" z
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ! Z1 B9 O) f6 q% k7 j\" A
   
8.         y=y+0.0009
   ; }+ D# I- I9 h! f9 P) Z# R- Q
9.     }, ' P* e: B' p# _# ~/ d4 P# d% H4 J
   
10.     x=x+0.0011 3 c. R8 k8 W# r1 \' |
    
11. }, , A. q! d9 }9 A/ |5 a, [# k: L
    
12. s;
    ' Z* ]9 @- a; M1 p8 b5 m% H1 ^
13. [sys::clock()-t]/1000;

复制代码

结果：
* t' _* V& Z! {" s! J1008606.64947441
0.734   //时间，秒
7 G4 r; ~' @! T5 ~. }
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
2 v% C, ?; |+ s; D
－－－－－－－
& K. ?- z' g- M8 z
# ^3 e8 y" F0 ~Forcal中还有一个函数sum专门进行这种计算：

1. mvar:& v# P\" @1 q. W1 J
   
2. t=sys::clock();
   5 Z- t/ o0 \) V; P+ o( Y8 E* h
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 \/ U! h, e' W# @; Y  y: [
4. sum["f",0,1,0.0011 : 1,2,0.0009];% n: }9 |6 d3 h4 w. y+ E, V# O\" w7 F
   
5. [sys::clock()-t]/1000;

复制代码

结果：
' z4 ?- [4 l1 @* Q2 Q5 x2 d( S1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
- H2 q2 l) G% O: q
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));6 G4 R; i! E+ m8 i! j, v+ A
   
2. tic
   \" r  y: y4 S3 R7 e& b, N+ M\" k
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   0 A\" [! w( r8 b' B- C2 r
4. sum(sum(arrayfun(f,x,y))). W9 Z* w$ M* ?
   
5. toc
   . [6 ~\" ?# w\" O( z2 r
6. 
   ) H' i1 M5 p/ _
7. ans =9 I3 ~+ o2 s+ l* S  m
   
8. 
   4 h1 g; H. C' h0 F, w% s
9.   1.0086e+006) J9 q9 k7 p\" I8 l
   
10. 
    ) ~8 J, \8 W0 f; O
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; _( `2 J/ l. E' f+ J5 I
4. mvar:( j( ^8 W! F$ D% e  G
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 1 G\\" c: p  @8 ]# G
8. t=clock(),. F4 J& G' I$ r4 f
   
10. oo{3 X: u' C$ }5 j' ]  `+ s
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],7 n, l) h# w5 ~; m
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. % a- }  ?2 p9 d2 C5 t8 p8 H
16. };$ u( m6 Z/ j& L5 f) I* n
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; _( `2 J/ l. E' f+ J5 I</span> mvar:<font class=\\"jammer\\">( j( ^8 W! F$ D% e  G</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\\">1 G\\" c: p  @8 ]# G</span> t=clock(),<font class=\\"jammer\\">. F4 J& G' I$ r4 f</font><br /> oo{<font class=\\"jammer\\">3 X: u' C$ }5 j' ]  `+ s</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">7 n, l) h# w5 ~; m</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">% a- }  ?2 p9 d2 C5 t8 p8 H</span> };<font class=\\"jammer\\">$ u( m6 Z/ j& L5 f) I* n</font><br /> [clock()-t]/1000;

结果：
0 n% r- r3 k" |- g- Y1008606.64947441
! N- ~+ r0 n: a3 d2 d9 l& x- A/ W0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
$ w0 e+ j5 P" [' {9 v) m6 h, [3 q" Q
9 e' w5 p( z0 C7 c. K, g" d5 c－－－－－－－－
4 ?- F( ?+ t0 N! ]1 {: e+ S2 o& y
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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