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

forcal

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

1. //用C++代码描述为：
   8 l0 w% O2 ]; N$ U: a! H6 z
2. s=0.0; 2 Z! u, I  R* c/ u2 P  A8 _, \& u+ a
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   # K4 X2 i8 E* P
4. {) w0 s( I. {\" ~) f% f) t
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)  `5 i# H) P3 i0 e' i1 @: W
   
6.    {
   & ?. @, G. A! m% ?9 `- J4 x& s
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( N3 r& ]7 `8 z0 m9 H4 c
8.    }
   ! f7 v5 t4 O4 X& d& ]
9. }

复制代码

Matlab代码：

1. tic
   & Y5 Q# i9 x* U% R+ P3 ?
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);. ]! x: ~; T/ k- `
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   # A! \0 O, k' h, c( {
4. toc- t$ [; {7 v+ z0 q- k1 R
   
5. 3 `& K' |6 c) h+ }! M
   
6. s =
   2 c6 S) b, n* V# C( p
7. 2 y9 D, S. j6 Q3 z
   
8.   1.0086e+006
   ( \, f6 l6 G; }- b: ^: [
9. 
   ' D* q' {2 h# K8 L9 V. l: ?1 N% N
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 3 x7 ~( |. V7 }* h
4. mvar:- Q' v+ f- k- W( P9 X9 _
   
6. t=clock(),
   
7. % K, ?  G\\" a* b
8. oo{
   
9. / R. a+ N\\" ~' s' v% H* B
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],9 V# t* O% X8 ^* x/ T/ j2 H: u
    
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]7 C) G  z4 C4 f' w  P
    
14. };, k# j  n2 L; |3 y' O9 i
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">3 x7 ~( |. V7 }* h</span> mvar:<font class=\\"jammer\\">- Q' v+ f- k- W( P9 X9 _</font><br /> t=clock(),<br /> <span style=\\"display:none\\">% K, ?  G\\" a* b</span> oo{<br /> <span style=\\"display:none\\">/ R. a+ N\\" ~' s' v% H* B</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">9 V# t* O% X8 ^* x/ T/ j2 H: u</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\\">7 C) G  z4 C4 f' w  P</font><br /> };<font class=\\"jammer\\">, k# j  n2 L; |3 y' O9 i</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
+ p' [3 F2 `' `. p9 v0.641
9 S  V+ D: U, l, d! `
Forcal比Matlab稍慢些。

8 I* R3 D4 h' \! k  j. o7 F% ?－－－－－－－－－－

, P6 K% x, T1 }1 |0 q  Z) _: C再看循环效率。
! F4 ^! |, Y. \/ A+ h5 h, }+ e
Matlab代码：

1. tic( \, I. `5 B8 x) g* O
   
2. s=0;
   * b; w9 v& C+ O8 ]$ w8 d- M
3. x = 0;1 R% @' u+ t- U# H0 u9 Q! V
   
4. y = 1;
   % I+ w) R: ]. _' P2 D! f1 w
5. while x<1   
   3 `3 Y\" e; J2 e5 q+ |
6.     while y<2;        
   3 Q- B7 R, K- E  l/ i8 r# A: k: f
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   7 f' Q\" Q+ b+ F/ i
8.         y = y+0.0009;        + J& s9 H  D4 w0 O  ]
   
9.     end
   ! c4 w& @+ h2 Z9 L( e4 s4 G
10.     x = x+0.0011;
    : j2 t. z4 i& ]  s, i; ~% p
11.     y = 1;
    ' k& J9 E* W+ J/ u
12. end* e4 d1 J8 c' l8 M: l! j; L0 N
    
13. s& ^9 t6 k9 P8 i\" Q  s
    
14. toc, R4 Z! q' G7 o) d0 n7 w+ D
    
15. $ q2 l- K. D1 X1 d
    
16. s =
    7 u7 q) W% a- X/ M
17. 0 R1 B* B& d3 |& M5 `6 a) h
    
18.   1.0086e+006  n  L+ w! g9 c+ H, s
    
19. / z. p6 p- {  e2 d7 S\" N4 O7 i
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   / K% L# `- x5 I2 [  y- F7 |
2. t=sys::clock();( q: }* O: o8 S* {- x9 c
   
3. s=0,x=0,
   , ]0 @1 @5 Z2 P, Y9 G! [% L1 m
4. while{x<=1, //while循环算法；
   . [) V( q) X\" g: `$ Y+ `
5.     y=1, ( V% D  W  O* H4 Q( r\" t
   
6.     while{y<=2, ; {0 S/ J\" U9 y+ x3 z+ C
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   , ~+ F% E9 {: V. ]- ]7 E
8.         y=y+0.0009
   5 @: c8 `, o7 s4 L  @& u
9.     }, . p+ p6 ~( B% @( _& L9 l
   
10.     x=x+0.0011
    ( _0 h: P5 v' F) |  c2 Z
11. },
    - k7 h: ?  L! N: }( o! ^& E
12. s;
    3 W8 h\" P% p6 Z0 E% |( F- ?
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
* P% m1 z  u& V  O2 ?6 A2 k8 W
7 J; ?9 h; n4 G5 d7 i+ E$ W# F+ }6 I我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
" k& W$ w$ H, n
$ w3 w# w# D1 {4 ?1 v. X－－－－－－－

Forcal中还有一个函数sum专门进行这种计算：

1. mvar:( q\" N% M6 F3 H1 w/ I0 ~
   
2. t=sys::clock();* G+ L7 f4 u: w) G/ V+ ?/ \0 N
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ( I$ g( y# b% d4 n9 D+ ?
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   1 t, v$ a. O3 G0 f& H0 ]
5. [sys::clock()-t]/1000;

复制代码

结果：
. G- a2 s# y, g5 X1008606.64947441
+ W: M' d- U& T9 X% T0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
5 s: N7 F0 ~- K) W# H3 k
" C$ d% ]9 g1 u- Nmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));6 ?' t4 }/ L$ E8 R
   
2. tic
   3 d  ?. v$ T- W\" m1 f4 t% n' D% W8 P
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);3 B; [& {: E% A3 R+ p+ W( O  @$ c/ F
   
4. sum(sum(arrayfun(f,x,y)))0 p3 d  E8 b2 o2 c
   
5. toc3 Z6 d# ?# S4 b/ t5 w( f% ]
   
6. : n. g' }+ S; A/ B0 K3 d6 c
   
7. ans =+ ]* D4 o$ [# o+ F6 a9 n# y' [
   
8. 
   + r+ f# x6 `+ E# z% d+ e6 z
9.   1.0086e+006( Y\" J) J; E2 c% b  M* O1 W! @
   
10. 
    + N+ f! p8 L8 J. p) D. h- v  L
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];- H; R( s, T\\" h. C
   
4. mvar:
   
5. , H! U* k# U! N% ~% r\\" q
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 8 r+ j, q& |+ Q& ^, i1 D9 P5 X
8. t=clock(),
   
9. 2 c$ [0 A% O5 w* l7 @( X
10. oo{
    
11. 8 i4 b3 |: i9 \3 S8 e
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],7 I, X* \. B/ q% t. u6 P
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 5 C6 }0 P! N3 K4 T2 l3 y
16. };5 `, g& a( `3 Y# T. n+ V
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">- H; R( s, T\\" h. C</font><br /> mvar:<br /> <span style=\\"display:none\\">, H! U* k# U! N% ~% r\\" q</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\\">8 r+ j, q& |+ Q& ^, i1 D9 P5 X</span> t=clock(),<br /> <span style=\\"display:none\\">2 c$ [0 A% O5 w* l7 @( X</span> oo{<br /> <span style=\\"display:none\\">8 i4 b3 |: i9 \3 S8 e</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">7 I, X* \. B/ q% t. u6 P</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">5 C6 }0 P! N3 K4 T2 l3 y</span> };<font class=\\"jammer\\">5 `, g& a( `3 Y# T. n+ V</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

8 d+ q$ f3 F5 D) b. c4 `& j+ y& `" p1 r－－－－－－－－
& _$ {/ u: |+ Z. Q# q+ e9 X6 G+ i- a
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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