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

forcal

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

1. //用C++代码描述为：
   3 Q6 p% r# U, C4 T5 Y
2. s=0.0;
   2 H0 L1 u4 S: y# J. s0 i! q8 E
3. for(x=0.0;x<=1.0;x=x+0.0011)
   - W. w\" b- T( A& s- h0 w/ c
4. {
     w9 V# f& ?( ~* D' W! H, A
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   $ E1 F- W6 b! l
6.    {2 u/ ^: W& h4 z% V# M3 `\" I5 n. E
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));# \; ~1 ]; i! [\" w
   
8.    }
   7 o% o* a; s& ~* w% [* |5 Q
9. }

复制代码

Matlab代码：

1. tic
   , Y0 u7 j. ?% w8 \2 S7 A\" q
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 @/ {# i: y4 A4 }7 M$ |4 U1 c
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))1 C' s2 N! ?5 T( n: {+ [2 b
   
4. toc# L- w, c7 D& I  h0 p8 d6 Y
   
5. 0 P\" h5 w) K1 {( ^& K\" U
   
6. s =; i$ W0 @% B7 W7 K# u( J
   
7. 
   8 [+ z% i\" e9 J( _
8.   1.0086e+006- a6 F; w! `, B3 F, u2 {0 w7 N\" ~7 f
   
9. 
   8 n! p+ d# h9 X+ [2 h& P
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. $ p1 b8 A/ |) q0 X2 |
4. mvar:
   
5. \\" N8 E+ W7 \+ g
6. t=clock(),
   
7.   D% D) H6 k\\" J9 N/ ]
8. oo{) O& y* o* X* D+ S7 l8 F4 b
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. ) l: k3 X* V' r) _6 J5 a: L% w
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. 9 ?% i9 H; k, I6 W
14. };
    
15. 3 |8 a\\" J' Y0 O: S; p9 I
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">$ p1 b8 A/ |) q0 X2 |</span> mvar:<br /> <span style=\\"display:none\\">\\" N8 E+ W7 \+ g</span> t=clock(),<br /> <span style=\\"display:none\\">  D% D) H6 k\\" J9 N/ ]</span> oo{<font class=\\"jammer\\">) O& y* o* X* D+ S7 l8 F4 b</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">) l: k3 X* V' r) _6 J5 a: L% w</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]<br /> <span style=\\"display:none\\">9 ?% i9 H; k, I6 W</span> };<br /> <span style=\\"display:none\\">3 |8 a\\" J' Y0 O: S; p9 I</span> [clock()-t]/1000;

结果：
1008606.64947441
8 n* \; }+ h' S) i+ H& F0.641
& ^! r( G: C4 {
  K# ]- D, |' ^1 V' i- y- }0 i5 NForcal比Matlab稍慢些。
% k0 J- g4 o. i! ?( U  b0 V
; f6 b2 @$ ^9 T" N/ Z－－－－－－－－－－
: G3 K$ |! n' |; S# r
再看循环效率。

& R/ c. q' o2 ^' d  _+ ?3 i, xMatlab代码：

1. tic
   ; E3 K* B0 }3 v# H4 w
2. s=0;
   # u3 x* U1 D; I0 s% a0 R$ {, V, _
3. x = 0;5 h0 j9 S  ^6 g  `& ~# q# g
   
4. y = 1;. V9 C: V8 y9 s+ V\" T8 L0 `
   
5. while x<1    ! ?; ]; a( C5 \1 t8 q) z9 d
   
6.     while y<2;        / a) f; J. Q0 }7 ?8 Q, t
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));' S$ r# u2 o, r1 Q8 p2 w3 Y$ S$ q
   
8.         y = y+0.0009;        
   0 h' N7 \6 x3 o6 h
9.     end6 \6 ~6 ^0 R7 a1 V
   
10.     x = x+0.0011;# S) S! {6 [4 ?# W( U, K\" }
    
11.     y = 1;/ \# V5 C9 J9 H0 d- z+ J' w3 W
    
12. end$ _+ X& X2 v  `' i# X# R, O$ ?% o
    
13. s
    # r' y: C\" m( p% {
14. toc
    & Q3 U' ~/ p/ _, A* q\" Q
15. 
    ( F) F6 K\" z\" c; F) B4 U
16. s =
    3 G0 l; P2 h$ i) q# `; K+ N
17. 
    % b2 M( p- Y# g. _
18.   1.0086e+006+ E6 `# H7 z3 n$ Z1 R& q' ]
    
19. 0 {4 W! `% q+ x+ b! E
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   * h# g7 U4 g; q& Y5 \7 T( d
2. t=sys::clock();
   \" }# b5 l4 X/ F) C; e3 O
3. s=0,x=0,
   % a! W' o/ X) q6 _6 k+ J
4. while{x<=1, //while循环算法； 3 j5 G! f* F) z+ M, e8 c# h2 Q
   
5.     y=1, . m& z9 @4 Z\" [- {  W1 D
   
6.     while{y<=2,
   4 n/ @7 X, w+ o) h1 Q
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   5 I  H: A4 ]+ S  S4 Y\" `
8.         y=y+0.0009 : f+ M+ V* G. K2 y) `# K
   
9.     },
   6 A- G, k3 ?- m7 s+ x
10.     x=x+0.0011
    : t; C; W\" p/ h6 W
11. }, ! ~2 {4 q, c: ?; q$ d& ~
    
12. s;\" Q% y3 k: l* z/ P
    
13. [sys::clock()-t]/1000;

复制代码

结果：
2 z) C, \0 T" O5 A8 b1008606.64947441
  y! t2 j0 i1 ~0.734   //时间，秒
0 ?0 {( U" v: L9 I9 W& M/ V# t
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

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

1. mvar:
   & [8 y; q; }1 P/ `. ^* K) ?
2. t=sys::clock();- c6 m6 Q& M3 U+ q8 A- I6 P  C
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 7 _2 g' D3 K9 Q0 o, f
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   6 {3 O+ J. O$ A$ q+ m$ h: t
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 J* d9 A\" j9 }# _
2. tic) ?. O; h) v\" J. N6 g
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ( a\" M, Q5 U* s) Y  |( l
4. sum(sum(arrayfun(f,x,y)))3 p+ x6 B2 s2 {( ]6 w0 a. V\" m$ S
   
5. toc
   4 [; ~2 B\" b3 i: {7 l
6. 
   , B, r- ~' r) x3 e+ Y6 z! j
7. ans =0 T' G5 D; c0 Y5 I  ?
   
8. 1 W1 S0 ^- V! P
   
9.   1.0086e+006
   - C0 ?, u' m: ]8 [2 ]
10. ; e8 w7 i; H3 ]- t! p5 l
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 2 i% `6 u1 b  @. R3 b
4. mvar:\\" }5 E/ n' o. n3 H\\" b- |  H
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. / @# d7 Q7 d0 I' V/ R* [* M
8. t=clock(),3 x( c9 p2 y4 {& j  ^( b4 [\\" P4 U
   
10. oo{
    
11. $ L0 _% q3 Q, U: ?% U
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. # L7 i& t( i  X+ d& {: F! a
14.     arrayfun[HFor("f"),x,y].Sum[]- R# T  U5 T$ t( ~. J
    
16. };
    
17. 5 T9 U3 F6 H+ A' G( c0 k. G
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">2 i% `6 u1 b  @. R3 b</span> mvar:<font class=\\"jammer\\">\\" }5 E/ n' o. n3 H\\" b- |  H</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\\">/ @# d7 Q7 d0 I' V/ R* [* M</span> t=clock(),<font class=\\"jammer\\">3 x( c9 p2 y4 {& j  ^( b4 [\\" P4 U</font><br /> oo{<br /> <span style=\\"display:none\\">$ L0 _% q3 Q, U: ?% U</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\"># L7 i& t( i  X+ d& {: F! a</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">- R# T  U5 T$ t( ~. J</font><br /> };<br /> <span style=\\"display:none\\">5 T9 U3 F6 H+ A' G( c0 k. G</span> [clock()-t]/1000;

结果：
" Q3 m, Q" k" q% v, k' z- X: \1 D1008606.64947441
; S  T0 e7 U# N- O, f, T( M0.735  秒
2 C2 Z. f4 e" }
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

' F  H' }; K7 l& o( }' x从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。