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

forcal

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

1. //用C++代码描述为：! N# a5 T; M2 _1 ?
   
2. s=0.0; ( g( _- A; H( w. {: |. c) o/ `4 @
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
     l* s) B: Y+ M# c6 K( g  Z1 D* V+ K# k
4. {+ [. w8 N: |, \+ H: e# g
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)$ t+ x* K) r/ ]  b$ _6 ?# _
   
6.    {
   ; r5 I. g4 B5 n\" O
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   $ `( W  m7 X7 F
8.    }
   6 m$ G! y+ X8 d# x; ]: I' B1 G
9. }

复制代码

Matlab代码：

1. tic4 j\" E- ~; H. E' m) P/ N
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   2 f8 G$ }# C* X( Z
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))) G& q1 `$ C8 N; c4 V, g2 r5 F
   
4. toc
   & V# P8 P* o8 _* s\" l0 b
5. ) e\" \+ c\" N1 f: M0 t1 k- B
   
6. s =
   5 J+ T4 ?\" b8 y0 U\" @4 D) W$ E
7. 
   ; Z# E2 f) d* K& \( c, g
8.   1.0086e+0066 m4 l9 e/ P- s% ^
   
9. 
   # G* b' e3 `% L9 ]. |9 ?
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 6 e% J# ^, }  y3 Z( q* p
4. mvar:- I& _8 K1 f- F
   
6. t=clock(),
   
7. , e# A\\" P+ ]5 l7 `8 z6 P3 J0 G
8. oo{
   
9. ; m' ^6 |& G; e6 \- A0 ^
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( |7 e( G' ^+ D, x2 M! J9 |
    
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]) M/ |3 ^. z0 n$ C4 J+ v: A
    
14. };
    
15. . A$ G% z- B% C' h9 D* y- v3 M
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">6 e% J# ^, }  y3 Z( q* p</span> mvar:<font class=\\"jammer\\">- I& _8 K1 f- F</font><br /> t=clock(),<br /> <span style=\\"display:none\\">, e# A\\" P+ ]5 l7 `8 z6 P3 J0 G</span> oo{<br /> <span style=\\"display:none\\">; m' ^6 |& G; e6 \- A0 ^</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( |7 e( G' ^+ D, x2 M! J9 |</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\\">) M/ |3 ^. z0 n$ C4 J+ v: A</font><br /> };<br /> <span style=\\"display:none\\">. A$ G% z- B% C' h9 D* y- v3 M</span> [clock()-t]/1000;

结果：
1008606.64947441
# c, w, B2 F8 S: j- f0.641

: U+ ?+ o% C" A+ l% a7 p7 yForcal比Matlab稍慢些。

+ \& A, ]2 D3 a* A2 b6 i－－－－－－－－－－

再看循环效率。
+ D9 [: ]0 @6 R, ]! s, C
Matlab代码：

1. tic
   ! r1 M6 k/ \  I% F0 g
2. s=0;7 b2 ~1 \2 O5 ], j1 @4 w8 ^- a
   
3. x = 0;
   6 P9 f+ i+ N% |! l. e2 g5 O& ?
4. y = 1;
   1 [) |$ Q; V, b/ g: A' n1 K/ K
5. while x<1    4 P( a& x3 }' e$ [( }1 }
   
6.     while y<2;        
   % r9 X+ B! q* {3 ^5 ~& T6 x
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 l0 V) ]+ Q; S+ {4 j% y
8.         y = y+0.0009;        # p( I  A4 N/ B2 k\" B- {
   
9.     end  H, G: u' E\" {( E; s- F+ e
   
10.     x = x+0.0011;
    4 [7 w6 x# T( b) \0 H) k* o& j
11.     y = 1;
    $ O$ w\" n! Y3 ]3 ]% [/ H! O
12. end
    & k2 n1 C8 f4 w7 G; K
13. s1 F7 C% s7 i; J: R+ V
    
14. toc; x7 b\" E- d: Q2 f: {
    
15. 
    ) O/ @\" g  w) L& b$ g& {
16. s =
    ' P9 E  B; R2 u6 g* N
17. 
    ( H5 U, c6 M: r, {, V. i
18.   1.0086e+006
    3 y7 U( W) _. Z  C1 S& {1 N2 w
19. . U8 X9 i! N3 ^! T) p
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   / e; ?1 a7 C$ [* n4 j. z
2. t=sys::clock();
   6 ^6 m5 P* q8 O\" z2 d6 V7 f
3. s=0,x=0,
   $ G0 d3 t$ f9 Z8 C' }\" O: G$ W6 r
4. while{x<=1, //while循环算法；
   # Q3 Y! {' K2 M1 }
5.     y=1,
   8 _1 P! u* n6 B; D$ t3 n
6.     while{y<=2, ! W7 }4 y  d2 {
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   8 N3 u; q1 |5 R) }
8.         y=y+0.0009 , a7 V9 s6 Y: d+ p
   
9.     }, 6 l9 o( F4 `1 ?5 l) t
   
10.     x=x+0.0011
    . v3 }$ C5 s1 U( k8 _
11. }, ! o! s; W% {4 a# E) x, U; q0 L+ G' Q
    
12. s;8 C, n( K) v# ^
    
13. [sys::clock()-t]/1000;

复制代码

结果：
/ l' y2 E5 H7 a1008606.64947441
  Z3 P5 ^2 B( v4 g( F0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

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

1. mvar:
   5 E4 J. l5 R1 @+ m
2. t=sys::clock();7 M  N+ s4 t7 I+ ?7 ?
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   + A+ x: W% V% A
4. sum["f",0,1,0.0011 : 1,2,0.0009];  l9 G\" E+ O0 y# E' j6 \
   
5. [sys::clock()-t]/1000;

复制代码

结果：
( P! O' Q3 j* s# U1008606.64947441
  t0 M! a7 b3 a" U: |0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
+ A& b4 k  E" T8 w. `4 }
. `# _0 N) j# @( p/ |$ o7 umatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 g8 b; F  I5 F1 y- G
2. tic
   # e. A5 [# g8 E; |  \4 `
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);$ H: ]4 X0 W3 ^2 R5 N' F* Q
   
4. sum(sum(arrayfun(f,x,y)))9 j8 v7 f4 J; W, x7 ^) a8 v
   
5. toc' o8 m% T$ C. O6 N+ Q& j\" N. Q
   
6. 
   6 c( p8 Y: j& ?- \& z7 a& o/ b4 N& p. Q
7. ans =
   & J, x4 r; w7 N1 a8 Q! T) k( C+ m' P
8. 
   0 o3 {( b7 w  g0 J* J
9.   1.0086e+006
   ! L0 U# V( ^* b
10. 
    ( V: y# x\" `6 c; I1 c' s* D7 _
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];; M6 g4 @! T+ \9 b) M
   
4. mvar:
   
5. ' \# g* h, O* u1 z* w7 @
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 9 m% V& X9 T5 E& q+ n3 B2 _
   
8. t=clock(),3 M3 w/ m' z; K; [/ W, {
   
10. oo{
    
11. 7 Q6 L3 ]! ~  X) x0 |
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],5 R7 v$ h' t1 j4 u4 d5 z7 {* |+ c
    
14.     arrayfun[HFor("f"),x,y].Sum[]\\" E+ E/ N% e5 o5 V; @
    
16. };; `4 h- x' }# [7 q. t% b
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">; M6 g4 @! T+ \9 b) M</font><br /> mvar:<br /> <span style=\\"display:none\\">' \# g* h, O* u1 z* w7 @</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">9 m% V& X9 T5 E& q+ n3 B2 _</font><br /> t=clock(),<font class=\\"jammer\\">3 M3 w/ m' z; K; [/ W, {</font><br /> oo{<br /> <span style=\\"display:none\\">7 Q6 L3 ]! ~  X) x0 |</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">5 R7 v$ h' t1 j4 u4 d5 z7 {* |+ c</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">\\" E+ E/ N% e5 o5 V; @</font><br /> };<font class=\\"jammer\\">; `4 h- x' }# [7 q. t% b</font><br /> [clock()-t]/1000;

结果：
9 x* }, x, W/ s0 L+ o1008606.64947441
0.735  秒
  X6 x4 ]6 L0 d4 c
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
( B; P: z1 O+ T, ^. z
－－－－－－－－

$ J! R3 F, c; @9 H) f从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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