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

forcal

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

1. //用C++代码描述为：
   3 D\" I& ~# @) |! `
2. s=0.0; . \: k. p' K( _3 s3 w  `& w6 x6 {
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   7 s7 I) v  e+ R
4. {\" k: O( g% P) o8 Q
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   $ R8 Z9 m# O; Z, f4 @
6.    {& W\" Z6 c\" ^/ j8 f  h
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   6 v) m9 y: r- m( `3 @
8.    }( B) X4 t\" I, t% o
   
9. }

复制代码

Matlab代码：

1. tic
   \" T! F& E! \. b$ c, f  U* X
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   . a( k- V% a2 o9 Z
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))( T3 W9 n6 P; B0 P/ o3 W  h1 Q. R  a+ ]
   
4. toc
   ( ]0 h\" r4 K1 p+ J
5. 
   # N\" K9 y% A+ P: A3 w' M; t
6. s =
   ; }/ ^$ ?7 ?( J% y
7. 
   ; U' C  W$ |$ ]' J/ {6 D
8.   1.0086e+006) _, g- ^  k& Y  D0 v' Q) @
   
9. . m9 H/ G2 m3 l
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 3 d$ B. U4 W8 G
4. mvar:* w/ K7 s  f; L. d. Z: D4 ?\\" t
   
6. t=clock(),
   
7. 9 c; H. Y8 Z4 Z3 z8 n
8. oo{
   
9. 1 i7 q/ K- x( p
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. ) H% l  v+ Q1 @# m- |: 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]7 }9 h& ^+ `; n! \  T3 \6 ^( U+ _
    
14. };) f& c6 w, [( N5 |. d% O- g
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">3 d$ B. U4 W8 G</span> mvar:<font class=\\"jammer\\">* w/ K7 s  f; L. d. Z: D4 ?\\" t</font><br /> t=clock(),<br /> <span style=\\"display:none\\">9 c; H. Y8 Z4 Z3 z8 n</span> oo{<br /> <span style=\\"display:none\\">1 i7 q/ K- x( p</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">) H% l  v+ Q1 @# m- |: a</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]<font class=\\"jammer\\">7 }9 h& ^+ `; n! \  T3 \6 ^( U+ _</font><br /> };<font class=\\"jammer\\">) f& c6 w, [( N5 |. d% O- g</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
1 `& Y, K" g0 z9 I  d; F, F1 }0.641
, ~/ g7 C2 Q, `, {8 V; {% t
Forcal比Matlab稍慢些。

－－－－－－－－－－
( N  p, [* U5 ?% A
6 I" O5 F6 r3 q9 c2 v+ ~" A再看循环效率。

, r' n7 z: u6 s4 W% R6 _Matlab代码：

1. tic
   ; {) F  T# D# }
2. s=0;  V) M2 S* V) h
   
3. x = 0;
   * |  L) B$ h4 u: F. j7 p9 ?- W9 b  d
4. y = 1;$ p8 }9 e9 o0 G2 m, a- T% z
   
5. while x<1   
   ) Z8 r/ z! w' [  M7 o
6.     while y<2;        * f4 A8 W. ?: h\" w. E6 t* S4 d/ N
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));; K3 h; P1 h0 z  {3 m6 z, Z5 ^
   
8.         y = y+0.0009;        
   1 s7 H\" J4 T2 @2 E% m7 @; k+ g
9.     end# I# K# V& P3 z' Y
   
10.     x = x+0.0011;  Z6 |8 K\" e1 o6 q: ^
    
11.     y = 1;
    8 i7 a7 Y' B; @* O# W7 M' g
12. end
    4 k. |) {\" x( k' t% x1 p
13. s
    . \\" t6 E$ q7 ?7 N\" r  L9 h3 y
14. toc
    & s& r\" m9 p/ e) E
15. 
    : l7 o4 l  ?! ^
16. s =
    - h8 Z/ U! m2 K6 D; W$ k
17. 
    8 j! }' ^3 y% r5 t
18.   1.0086e+006
      B+ Z: h: q1 {4 @
19. - _! b; s& k% I, ^! g
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:! V\" V6 |6 w  a4 v2 k+ ~
   
2. t=sys::clock();8 Z8 Q3 H\" A0 f
   
3. s=0,x=0,
   0 ]& a$ t! U& \, ?
4. while{x<=1, //while循环算法； % |6 R7 i8 M$ w( n6 D$ Z! C# }
   
5.     y=1, + n1 ?0 e2 i) D- e# F0 ~
   
6.     while{y<=2,
   6 }$ J* B1 j7 ]- K  d
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), : G: `1 X' P6 I4 o\" |8 o6 O
   
8.         y=y+0.0009
   . T+ {, ]% P5 I8 p1 x; b8 [: ~
9.     },
   ; F2 b7 g+ H4 A+ P\" F1 x
10.     x=x+0.0011 ! Z# v8 k, w# J# h
    
11. },
    * n9 g\" Z, J7 z
12. s;
    4 s/ _8 w& J! H: i
13. [sys::clock()-t]/1000;

复制代码

结果：
$ A2 p! m% H  w- P. v: |- y9 e" d1008606.64947441
8 T- R, l# B" B0 }7 }# R6 g% \* T' @0.734   //时间，秒
. T' Z& F+ F% p& _8 O8 a8 p
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
% D2 Y! K' M1 k$ g" G
7 B7 V- s9 Q+ g% g' g" ?4 J－－－－－－－
: Z6 i1 {: W; W* n
( w8 U. y& }. `2 l' r4 ^Forcal中还有一个函数sum专门进行这种计算：

1. mvar:( z  I; D+ d& [  v  U# s
   
2. t=sys::clock();) h' }9 O- z8 C- {9 Y3 T
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); * [. ?+ o* D& C
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];( N7 _( |! ?& ]* s. b+ S
   
5. [sys::clock()-t]/1000;

复制代码

结果：
& Z: Q6 h" K8 K$ p8 D  L" U- p1008606.64947441
2 y0 N: j- p0 f1 Y6 G7 |0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
2 _. ?6 i# e: F
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));( r  z8 T5 `0 \7 q+ B' t1 d. V\" v
   
2. tic
   2 r2 n7 T6 g, I* L) W9 o
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   1 f& Q$ p$ w5 j
4. sum(sum(arrayfun(f,x,y)))) }5 d\" Z6 K7 d& j% ~5 |4 D
   
5. toc' C( V$ a6 N0 ]5 E+ w
   
6. 
   & G8 {  @- D+ m
7. ans =5 x# ?6 i  r' ?0 J2 k
   
8. 0 @! m' A* K6 I
   
9.   1.0086e+006
   , c; w' ^9 f. Q& y
10. 
      q7 ?9 o, G& R' j7 L2 x& G
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];8 d$ a. Z7 R  C$ X5 r  ]# N
   
4. mvar:
   
5. ; G* u5 `. u; W- W( T
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ; J0 s$ C/ r% o# ?
   
8. t=clock(),
   
9. - ^1 V! b! P# A6 L) {( V2 b2 p
10. oo{
    
11. * F% y8 a7 T) |- C
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 1 k( N' M4 i8 L' Y
14.     arrayfun[HFor("f"),x,y].Sum[]: L3 N8 W; {$ Y' t( l
    
16. };% g/ ~& J' J# l
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">8 d$ a. Z7 R  C$ X5 r  ]# N</font><br /> mvar:<br /> <span style=\\"display:none\\">; G* u5 `. u; W- W( T</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">; J0 s$ C/ r% o# ?</font><br /> t=clock(),<br /> <span style=\\"display:none\\">- ^1 V! b! P# A6 L) {( V2 b2 p</span> oo{<br /> <span style=\\"display:none\\">* F% y8 a7 T) |- C</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">1 k( N' M4 i8 L' Y</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">: L3 N8 W; {$ Y' t( l</font><br /> };<font class=\\"jammer\\">% g/ ~& J' J# l</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
$ y6 L7 W! l& A9 Y0.735  秒

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

－－－－－－－－
5 |, x( ~$ ?" ^
2 E! u* W. K" e从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。