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

forcal

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

1. //用C++代码描述为：
   5 Z3 Y  t3 a8 l& p3 p7 U
2. s=0.0; 0 I! W# s% w+ k
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   , L4 U1 W; y\" @5 q6 R1 E3 l  Z
4. {\" I5 {% P1 k' z$ T# S8 g
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   4 r) D( s, ^* @
6.    {3 n0 G2 j; V0 u4 c- Q1 y
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & U2 I- d. i' y: x0 y
8.    }
   8 Z- h. T, y) E8 ?1 u+ h, n
9. }

复制代码

Matlab代码：

1. tic
   ( J' [+ R! z9 n# r! \
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);7 v: q9 H* ^& h\" I0 `1 m% ~- z8 _
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))3 \  D9 \4 _$ M+ L
   
4. toc
   1 K; {# B) u\" z9 L+ O! f
5. 
   4 P2 z/ k! v- U# o/ e% `  s
6. s =' [\" v6 X, w# T* T+ u( G* z; l
   
7. 2 k) ]  k& n6 S
   
8.   1.0086e+0066 s+ d6 I8 J$ W! J
   
9. 7 a/ k, \5 j- P: O( B
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];% }8 b2 U& @8 d. D
   
4. mvar:% X3 a- U' \2 F2 p# x
   
6. t=clock(),
   
7. , b' W: A9 I+ H) p* j
8. oo{4 R& W2 B\\" G* K! d
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. ' c' m  v+ `: H1 T8 J
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.   _. _1 h  i7 c: D$ j9 T
14. };( z0 J: f6 x/ @; R
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">% }8 b2 U& @8 d. D</font><br /> mvar:<font class=\\"jammer\\">% X3 a- U' \2 F2 p# x</font><br /> t=clock(),<br /> <span style=\\"display:none\\">, b' W: A9 I+ H) p* j</span> oo{<font class=\\"jammer\\">4 R& W2 B\\" G* K! d</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">' c' m  v+ `: H1 T8 J</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\\">  _. _1 h  i7 c: D$ j9 T</span> };<font class=\\"jammer\\">( z0 J: f6 x/ @; R</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
& U. |9 H; B* p" f2 b0.641
3 j6 H/ q) G% t5 T1 T
! `4 Z8 T3 k" Z; S8 S; JForcal比Matlab稍慢些。

$ M' W- y5 ~$ g" [－－－－－－－－－－
3 H& u* K" u+ X3 I# l. o9 e9 R+ H
再看循环效率。
8 e% y  f5 C0 \1 i/ R
, f' R; f# H( u4 _Matlab代码：

1. tic/ C4 A- _9 G* e+ G; I0 y  y
   
2. s=0;\" P4 q\" k6 K7 M5 B2 }( @! c  R# I
   
3. x = 0;$ M7 M! J) |5 J$ E0 A( {+ O6 X
   
4. y = 1;
   % k  M! {4 V1 e$ A; J4 |( q) v- E9 \
5. while x<1    8 l' M7 S0 T& g\" S4 m1 n
   
6.     while y<2;        
   % o# a- _7 u; B\" K1 n3 L' R
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));& f$ k% h; y9 W4 r! r
   
8.         y = y+0.0009;        ' F, o1 s& Q  z& d\" N+ q. x
   
9.     end
   8 B) f\" p3 g& B: t$ ^
10.     x = x+0.0011;# A6 F; K& H7 S# r$ ?4 s/ m% T7 W7 O
    
11.     y = 1;' X. `$ ^* a\" ^8 k5 _
    
12. end
    / z, I+ m0 R2 @% q4 Q, `# p
13. s
    ( L& J/ s- E2 B5 D
14. toc% s& `& J8 w' f; `0 P8 B- Z
    
15. 0 h. z& @7 n* ~
    
16. s =3 c# a: M: o' |3 k' {
    
17. 
    + S& R& A9 S% _7 m
18.   1.0086e+006% g# I2 K+ W5 S# {3 F; U8 Z. C
    
19. 
    5 u& K4 Y5 m! s; n
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:% y( M  ~8 _# {
   
2. t=sys::clock();
   / U& C2 t- d+ X' w  A4 n( W, w
3. s=0,x=0, 8 D+ G% N0 {! y- K6 x6 O3 M
   
4. while{x<=1, //while循环算法；
     P  ]; ^7 b2 ~% I9 q
5.     y=1, 9 H* Y2 F/ Z% l) {& ]/ Y% m' ?- v
   
6.     while{y<=2,   N; s! N! k4 }$ y& h, v
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   & t' l\" W  }9 e( V. {  [\" {9 C
8.         y=y+0.0009 ' K/ w6 s- s* w: A3 q
   
9.     }, 3 B4 ?( {4 A+ ]
   
10.     x=x+0.0011
    7 Z; w- v# S9 O# i
11. },
    ( e# [% [9 Q! p8 R$ y7 D3 O, s
12. s;
    . A, B3 @, T* v0 i5 y/ m
13. [sys::clock()-t]/1000;

复制代码

结果：
5 v4 w. |+ J2 {; f* {1008606.64947441
( D1 g% _2 j8 t& A$ H0.734   //时间，秒

; a7 s6 m5 K% U( L) f( u  X% E5 b我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

5 A0 ^2 B6 T. q& O－－－－－－－

3 w: O" P% q3 Z) N* S* mForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   4 t. w6 I& b' s
2. t=sys::clock();. J  ^* n# i3 R/ _. e8 r
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 h& G) d- C! x' i! b0 O
4. sum["f",0,1,0.0011 : 1,2,0.0009];. ~5 v3 W. N( X
   
5. [sys::clock()-t]/1000;

复制代码

结果：
6 O, P( v: f! S/ m* }8 T& u# M1008606.64947441
0.719   //时间，秒

forcal

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

' F; e/ e$ G0 J# o3 k: gmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));1 R( z\" t) S( {
   
2. tic
     c6 X4 F0 m% _
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);8 o$ t+ D( L( a. \  ]! J9 P0 ?6 O
   
4. sum(sum(arrayfun(f,x,y))), B% U9 K3 G\" v) s\" X. X3 m( h# }5 s4 O
   
5. toc# f9 |; x  Y0 P
   
6. 
   . K3 D. M/ t( n
7. ans =( x3 c6 q8 o) p& G/ @4 y1 E* N\" J
   
8. 9 f+ G$ |6 }) h
   
9.   1.0086e+006
   & H: |$ R; a; p7 E
10. 7 c  J, E' y6 e& d3 Q
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];7 W+ s\\" x3 y! f/ \* w
   
4. mvar:% j) R/ N& X4 m3 Z, b
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 2 a$ V' Q% D9 Z% ]
8. t=clock(),8 {% {% S) o\\" E  U3 P$ _  p. M
   
10. oo{
    
11. . b, c+ Y) j0 T; d7 S
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],2 I$ u; `. B! c0 I, r
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. . r- ?6 d% X. B: r, _* w) b7 b. t
16. };
    
17. * J+ v; v* E0 [* Y
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">7 W+ s\\" x3 y! f/ \* w</font><br /> mvar:<font class=\\"jammer\\">% j) R/ N& X4 m3 Z, b</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\\">2 a$ V' Q% D9 Z% ]</span> t=clock(),<font class=\\"jammer\\">8 {% {% S) o\\" E  U3 P$ _  p. M</font><br /> oo{<br /> <span style=\\"display:none\\">. b, c+ Y) j0 T; d7 S</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">2 I$ u; `. B! c0 I, r</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">. r- ?6 d% X. B: r, _* w) b7 b. t</span> };<br /> <span style=\\"display:none\\">* J+ v; v* E0 [* Y</span> [clock()-t]/1000;

结果：
1008606.64947441
% P+ A: p' R: z5 {; `0.735  秒

; `, v3 t- y) Z* J1 v可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
) H" W. N2 ], v$ I$ |# ?& S$ A8 ?+ D
－－－－－－－－
- J$ D' [$ E! `3 j% c
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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