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

forcal

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

1. //用C++代码描述为：
   7 e8 E1 W) G8 G/ H  A
2. s=0.0; 9 r& o# D! b& E4 V0 J; L0 L
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   ! _) ?- R/ X  o. E8 L5 ^+ x
4. {( L. f! O, E9 I4 l- Z) L+ U. T8 `
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   * }+ `\" A7 a0 D2 |6 ^5 W
6.    {
   9 q  J, X( B- C6 h0 j  R2 g, L
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));( F2 S: _' y1 C9 q& t% u
   
8.    }
   9 G7 u5 c7 p/ W5 A; q; Y
9. }

复制代码

Matlab代码：

1. tic
   ) n; ?% r8 v* |4 f+ }8 f0 a
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);$ q\" m! [& [9 u8 z\" P
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   / X/ M* Q1 C( N: ^9 Y
4. toc) O) ?+ @* {9 z/ f9 D
   
5. & X: b& S- T- l& @
   
6. s =  K6 Q& v, L/ C- q
   
7. & T( S2 C3 S; B3 ?' ]. f4 U% V# q9 v
   
8.   1.0086e+006* Z  f\" a( j* z' K# A& y2 K
   
9. 
   % m& E- e/ G. }+ u
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];8 S. g7 v3 C6 S$ c$ Q% P8 r
   
4. mvar:
   
5. 7 e4 P1 V* E; I' B$ K# \0 [# s
6. t=clock(),) k) _% Z9 l0 B! m. |6 S
   
8. oo{! h% k% u' s  ~4 k% y& m
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 9 ?% q\\" h4 R7 @* y
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. . E8 R8 b4 z6 }) V
14. };* G; Z# Q6 V- h1 S$ U
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">8 S. g7 v3 C6 S$ c$ Q% P8 r</font><br /> mvar:<br /> <span style=\\"display:none\\">7 e4 P1 V* E; I' B$ K# \0 [# s</span> t=clock(),<font class=\\"jammer\\">) k) _% Z9 l0 B! m. |6 S</font><br /> oo{<font class=\\"jammer\\">! h% k% u' s  ~4 k% y& m</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">9 ?% q\\" h4 R7 @* y</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\\">. E8 R8 b4 z6 }) V</span> };<font class=\\"jammer\\">* G; Z# Q6 V- h1 S$ U</font><br /> [clock()-t]/1000;

结果：
: h! l8 j4 P  p, o, w% \1008606.64947441
6 K, @0 W; ]3 U1 R% f6 s8 u0.641
2 f7 v- g1 S- ~
Forcal比Matlab稍慢些。

+ I2 n% b. j: Z! B' i$ v4 s－－－－－－－－－－
; Q( y+ @3 P& u0 [
: B3 S7 F$ D4 J" v1 B再看循环效率。

; c" c- Y& ~% [6 T8 V6 b6 QMatlab代码：

1. tic
   3 p  L* a2 W7 r6 p\" e0 h) _; W, O6 C
2. s=0;4 }+ h7 O. p1 V
   
3. x = 0;9 O+ b: F; W; X
   
4. y = 1;
   ! ]$ |2 s5 c4 c. p( a8 f
5. while x<1   
   + ~/ R5 |; u. U# y6 V
6.     while y<2;        6 \( T% Q5 n) y5 p0 @: O
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; ~0 |8 Y1 z8 {  F
8.         y = y+0.0009;        : q\" Q7 B3 q3 j& E9 Z: K4 |
   
9.     end
   ( L6 M. d0 ]2 O% c+ Y. Z
10.     x = x+0.0011;5 f% ?. m$ k& |$ o
    
11.     y = 1;
    $ u' _/ V8 F6 v% ], E3 {
12. end
    6 q. e$ D! x* W2 k
13. s
    2 J+ O! B- N4 V5 ~) }) S& _
14. toc5 W8 @1 }7 ?* u6 n; r9 p: V; {
    
15. 
    3 k7 n  s; e! a7 u8 k& F
16. s =
    5 i8 J5 O4 l7 k6 X( w) N+ N9 ]
17. & t2 d+ @' j/ |1 U& N
    
18.   1.0086e+0061 D+ u. B$ [5 d! {
    
19. 
    % Q7 x) S: Z  z4 n+ N( G7 m
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:* v# K! d2 t7 P+ O
   
2. t=sys::clock();: H0 `0 m' D# J. y  @3 A) ^( w
   
3. s=0,x=0, 2 c  g7 f: ~$ ^, W5 B\" ?
   
4. while{x<=1, //while循环算法； ' @/ F- F5 v( T+ W: S  E
   
5.     y=1, + H1 c2 {4 D0 |( {9 X4 d8 D
   
6.     while{y<=2, ( Q  {\" U* y' L5 \
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), $ _& @# V$ n( j2 ?
   
8.         y=y+0.0009   a4 `1 [8 u( I1 N, _$ v( ~
   
9.     },
   ; D* o/ _  `* Q' _# A1 Q
10.     x=x+0.0011 , k% H  P/ ]) S! p# C
    
11. },
    8 v5 q% L( B+ f
12. s;8 u* r6 {2 L, K6 B\" J4 A! h4 g
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
- h/ y) K* u, O# ~0 g0.734   //时间，秒

8 j6 ]4 S7 Y& _' U% Y3 X. G2 v我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

+ m7 V1 g7 s+ p$ o* L: i* |1 Z4 B－－－－－－－
$ s0 E6 H/ }) X' {. @7 E. Y- c
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:( v( W* K8 ^6 p4 L\" y
   
2. t=sys::clock();$ g% _1 v$ j0 k/ m; t
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 t. h% ?\" p: T7 \
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   5 a* |2 D$ O$ n& s( H/ h! l' M
5. [sys::clock()-t]/1000;

复制代码

结果：
1 f8 z9 @* h9 d1 a5 g6 G1008606.64947441
4 h/ v( g: {5 z4 w) x0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % y6 b( X' v: g  l) C
2. tic  r8 {4 V4 T+ i5 m\" f
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * [  z% n3 T9 X* m
4. sum(sum(arrayfun(f,x,y)))
   ( t7 |7 X$ a9 x7 o7 H2 |9 W$ }5 j0 {9 _
5. toc
   4 U, C  k( E8 A
6. & q% ^9 a$ v; J
   
7. ans =8 c2 f  r1 P. E
   
8.   |0 t4 @6 `7 ~9 K
   
9.   1.0086e+006
   - K+ a# `( L! @' w5 ?- O: U
10. : O1 ~% |# W6 N, Y
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];' I1 t* W3 d- V3 N2 y
   
4. mvar:/ d5 Q0 Q3 w+ L6 z
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); \\" C\\" z1 ]2 M8 E0 J\\" V
   
8. t=clock(),
   
9. , M' Y7 _2 N8 q6 W% k% u
10. oo{6 a* A4 M( ?2 |; `. W
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. . P. W  B& K2 I- L( I3 G
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 1 T! R- C1 B/ k8 d, Z) n
16. };
    
17. 2 a: U. w, n) [& d0 p* \
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">' I1 t* W3 d- V3 N2 y</font><br /> mvar:<font class=\\"jammer\\">/ d5 Q0 Q3 w+ L6 z</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">\\" C\\" z1 ]2 M8 E0 J\\" V</font><br /> t=clock(),<br /> <span style=\\"display:none\\">, M' Y7 _2 N8 q6 W% k% u</span> oo{<font class=\\"jammer\\">6 a* A4 M( ?2 |; `. W</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">. P. W  B& K2 I- L( I3 G</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">1 T! R- C1 B/ k8 d, Z) n</span> };<br /> <span style=\\"display:none\\">2 a: U. w, n) [& d0 p* \</span> [clock()-t]/1000;

结果：
1008606.64947441
6 ?; l- t: X% T' q( j8 {0.735  秒
9 ~6 S# i( ?/ M5 P: E+ E* v; w, ]. V
5 }# [2 S$ f# k7 ?" M% o可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
; [+ j3 u3 w: }0 R4 M
& {' h8 |  x! E+ j: x$ f8 _－－－－－－－－
/ r; R4 w! w' ?3 l
8 [( _* F6 [4 x/ I0 M0 Y从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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