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

forcal

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

1. //用C++代码描述为：) j4 U\" ^5 k7 ?$ ~6 U
   
2. s=0.0;
   + c$ }+ V; c- @8 g
3. for(x=0.0;x<=1.0;x=x+0.0011)
   1 X2 w+ I2 ^. t6 R& L4 o/ }
4. {0 d3 G. r5 T6 d, n
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   3 T9 @* k, p  R1 [6 A; H
6.    {
   \" G. Z- |$ ]& t. x6 T0 z+ Z
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 ^* r! A# \( t$ b
8.    }7 H- J9 I* K- h7 T+ B9 e
   
9. }

复制代码

Matlab代码：

1. tic
   0 W' y. f7 R/ B6 }
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * L/ T* n0 g; @( M
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   8 O0 b3 m0 s$ v4 Q/ a$ E; T
4. toc
   , y4 u* g; P8 B9 k, o! L: }* p
5. $ Y6 {/ C  u+ I  |% H9 d) T
   
6. s =1 e6 H! n* d3 i5 ?7 o+ {) H1 F8 _
   
7. 
   ; {3 J! D5 S: K& u$ y3 v
8.   1.0086e+006$ |. I5 K( s7 [7 v1 k
   
9. ( y4 ?: B2 t+ [1 i
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 5 z& J: j, M5 O* x2 t; b$ X$ c* h
4. mvar:, z6 R2 h0 M' W& @  h) v* s% P
   
6. t=clock(),0 U: h( @; Q# m( c
   
8. oo{% M' d\\" ]% [4 F! @. ?
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( H4 q9 k( G# 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]% i3 K1 k; L  c) b3 d1 I
    
14. };0 f: B2 k4 C; o9 }  K8 V% Z( D: F
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">5 z& J: j, M5 O* x2 t; b$ X$ c* h</span> mvar:<font class=\\"jammer\\">, z6 R2 h0 M' W& @  h) v* s% P</font><br /> t=clock(),<font class=\\"jammer\\">0 U: h( @; Q# m( c</font><br /> oo{<font class=\\"jammer\\">% M' d\\" ]% [4 F! @. ?</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( H4 q9 k( G# Y</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\\">% i3 K1 k; L  c) b3 d1 I</font><br /> };<font class=\\"jammer\\">0 f: B2 k4 C; o9 }  K8 V% Z( D: F</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
* F4 c+ o; o# j) t$ u- D4 t0.641

) m+ p3 h: Z9 p  C3 oForcal比Matlab稍慢些。

' o* a; s* |8 G) `－－－－－－－－－－

再看循环效率。

' m" {  ?8 c0 M/ FMatlab代码：

1. tic' b9 {( \3 @1 }9 S0 w
   
2. s=0;
   ; w3 }; ~3 P4 B! ]( ^& g
3. x = 0;
   . _1 ~; C( \% O\" w2 q
4. y = 1;
   4 S/ B& o  I% K  o* W& n* c* G- J
5. while x<1   
   * y- K% [5 Z/ |3 g; |( g
6.     while y<2;        
   ! b2 ]0 U3 X) n  f
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 p# l3 @2 V( P8 ^# |4 e: b  a; t
8.         y = y+0.0009;        8 |8 b\" {, \  P: d
   
9.     end
   - S; U2 o4 F' Q  T5 S& K. V) H5 t
10.     x = x+0.0011;
    8 w- S- o! J9 w6 i; z4 t
11.     y = 1;
    0 ^) T\" e( w% Y: J- J, ^
12. end$ B0 b$ o6 W* A; a* y6 R+ e
    
13. s
    * [* i+ x( P* S3 O& S
14. toc
    * Y! \2 f2 r# L8 d% y
15. : P8 K$ d# M! U' U5 Z( I. N
    
16. s =3 D6 u\" T; a' o# p) V0 g  H& R
    
17. 
    - N3 [' I7 E1 E- _7 q
18.   1.0086e+006$ r8 z# u0 x0 b9 }1 i) @, J
    
19. 
    4 p# p/ ^3 v8 J- W/ {
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   # @3 ~  d3 x$ ?4 I9 G\" L' G+ q) k
2. t=sys::clock();+ X6 ^& i7 ^; a: m- R! S- y
   
3. s=0,x=0,
   # m! k/ o+ {/ w\" ^/ ^3 p: X  G
4. while{x<=1, //while循环算法； 9 Z2 P7 P) n$ O, m4 f! r1 Q
   
5.     y=1, ! c, o9 W( z8 W4 Q! R3 \/ c) W
   
6.     while{y<=2, ) ?8 t( \# K1 X
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), & U7 Y  t3 y* M8 M! Y\" t+ C
   
8.         y=y+0.0009 ( M1 W6 `7 K2 w5 G
   
9.     },
   8 v/ {2 v; q' g0 @7 ^
10.     x=x+0.0011 7 n\" q$ o0 T8 P' B) [4 X
    
11. }, . v\" j; ]# a6 j
    
12. s;
    . j3 o! c. ^8 K1 o$ N& o
13. [sys::clock()-t]/1000;

复制代码

结果：
( M& ]' M  I8 s4 K# u5 Q1008606.64947441
0.734   //时间，秒

! ^  F/ ?3 w8 e5 V我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
$ J+ n" J8 H$ \  |' Z" j- N
/ Z0 k2 T: S4 _0 [, ~－－－－－－－
  D9 T! }. x) s7 W9 W
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   & {. U/ ~6 @7 `8 n4 y4 [
2. t=sys::clock();
   % r6 C* `6 U/ c- d/ y9 z
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! [1 I\" ~+ ~, q. {6 w3 \. r0 T& p7 Z
4. sum["f",0,1,0.0011 : 1,2,0.0009];8 S$ O0 E0 I+ Q4 p9 k# u
   
5. [sys::clock()-t]/1000;

复制代码

结果：
/ H- D# b+ o) S1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
) f8 A4 m! [' T6 S
$ W2 [7 l- i: e& a) Y. A; Omatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));. _  X1 C: J$ p( N8 g+ W5 O& J\" M0 J
   
2. tic
   ' z% C% X+ X/ V. r2 x4 C/ g
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, \- T, {5 I3 t: L
   
4. sum(sum(arrayfun(f,x,y)))
   ) c* ^# M( R) R\" j! c
5. toc2 V. F( F: e3 x! R6 O) ?
   
6. 
   $ ]3 H; H; e3 R3 p
7. ans =
   - r# M; `! h6 D  }/ R
8. 
   / ]9 A9 z5 w4 u6 W5 V( [9 n) F
9.   1.0086e+006& G2 D3 @; i7 Y
   
10. 
    ; a& V' G( `\" ~* G: f
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 8 }8 x! T4 A\\" f: k0 M0 ?1 O' R0 q/ [
4. mvar:6 X2 c  b+ h* n  A+ [+ L# m
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 3 d4 L  t1 Q: Z2 G$ n7 M5 |
   
8. t=clock(),
   
9. 6 @/ B4 B1 A6 m  S2 u
10. oo{, h\\" t& T+ c6 |) Q* g! W+ G
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. \\" w, @/ ~: C) H  `) ~
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. : b; T9 ]% @! H7 ~
16. };3 [6 X1 d0 v: n; c
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">8 }8 x! T4 A\\" f: k0 M0 ?1 O' R0 q/ [</span> mvar:<font class=\\"jammer\\">6 X2 c  b+ h* n  A+ [+ L# m</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">3 d4 L  t1 Q: Z2 G$ n7 M5 |</font><br /> t=clock(),<br /> <span style=\\"display:none\\">6 @/ B4 B1 A6 m  S2 u</span> oo{<font class=\\"jammer\\">, h\\" t& T+ c6 |) Q* g! W+ G</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">\\" w, @/ ~: C) H  `) ~</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">: b; T9 ]% @! H7 ~</span> };<font class=\\"jammer\\">3 [6 X1 d0 v: n; c</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
; \" ]& W3 k4 c6 |- i0.735  秒

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

－－－－－－－－
3 D9 k' _$ b) G- y* M
7 g' o" `' ?# C; C6 T从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。