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

forcal

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

1. //用C++代码描述为：$ S# K. F/ Z* ~1 E- z
   
2. s=0.0; 1 D( S6 X$ x) c
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   5 F& m! c: B* Y
4. {
   ' N( l$ i! r& `0 s
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   5 W4 n/ M$ u% [6 J! o, A
6.    {
   & g$ Y1 V  D, D/ ?( d$ n+ `
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! d4 o3 u. |; D9 n! E: w. P* J6 i) I& y3 u
   
8.    }
   6 S* m# ?) z  k. a
9. }

复制代码

Matlab代码：

1. tic
   : [( h9 R. t' ~0 I
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   5 y# l8 Z. R& e. Q$ n
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))6 r, j/ N, l. M+ z+ Q
   
4. toc
   / F2 S; p/ h$ z- X1 Y3 r
5. 0 d6 ]  }. c% Z1 S9 _$ g
   
6. s =
   % s- U: B- T, [
7.   i  }, ^! i- i% j8 v
   
8.   1.0086e+006: l& x/ @5 Z4 g& l* R* V
   
9. 
   ; S9 O. S8 R& @3 G9 }( F) \
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. : `/ A( _$ a- K' ^, O2 E! |4 P
4. mvar:
   
5. ; U+ Z. u5 t) p! q\\" }
6. t=clock(),
   
7. & @% o* w& v1 M$ I, q/ L4 [* D
8. oo{
   
9. ( t: m# u, l( e) {: [
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],8 s8 o1 |) M4 `
    
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. 7 a: O* T, J9 N+ a, c) W9 {$ L! u
14. };
    
15. 9 i/ j. |) H) ^( o- y3 N8 r
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">: `/ A( _$ a- K' ^, O2 E! |4 P</span> mvar:<br /> <span style=\\"display:none\\">; U+ Z. u5 t) p! q\\" }</span> t=clock(),<br /> <span style=\\"display:none\\">& @% o* w& v1 M$ I, q/ L4 [* D</span> oo{<br /> <span style=\\"display:none\\">( t: m# u, l( e) {: [</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">8 s8 o1 |) M4 `</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]<br /> <span style=\\"display:none\\">7 a: O* T, J9 N+ a, c) W9 {$ L! u</span> };<br /> <span style=\\"display:none\\">9 i/ j. |) H) ^( o- y3 N8 r</span> [clock()-t]/1000;

结果：
  k. [" a* Y3 H9 L. t  Q3 u1008606.64947441
+ G. e9 S! r0 \0.641

2 X" `$ E1 I2 T+ W. B+ M1 ^7 |1 x5 oForcal比Matlab稍慢些。

' J( F) c" Z$ o$ O! h% @9 _－－－－－－－－－－

& Z! [8 R" D% Q8 \% J, U+ P再看循环效率。
6 o: E/ ?' e5 b( O8 ^( q
Matlab代码：

1. tic7 N3 ]& j  t! F
   
2. s=0;9 N4 i4 R( }0 z, M8 J
   
3. x = 0;
   4 S, Q6 I# G* Q! C8 o& ]\" {! U
4. y = 1;
   - s% L5 ]$ l, w$ s0 q2 W: v; E
5. while x<1   
   4 @2 c3 l* O7 T% c$ a/ j: L
6.     while y<2;        % G, `% o/ `/ i\" u  G( }, |. a9 J
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));: r1 p6 m. V' n1 b
   
8.         y = y+0.0009;        5 ~( }- ]% t/ {% L0 V$ |$ z( c7 P
   
9.     end
   - L\" v% F, P7 r6 ^0 k$ ^0 ^4 q
10.     x = x+0.0011;9 Z, j- w- }+ |
    
11.     y = 1;7 V3 U5 p7 e9 I% x) j- U
    
12. end1 I\" c3 U6 H& G1 g/ E0 e7 N
    
13. s
    4 B( ~2 ^! q% F, `) O. V% r
14. toc
    2 N3 H0 S& b1 y
15. 1 Y  {- h6 A& w+ Q
    
16. s =% r& n9 V+ r8 E. m
    
17. 
    % i  M  L8 A2 `6 j: g( L: k0 ^7 J
18.   1.0086e+006
    7 s4 Q\" N8 X. C8 l
19. % T, k6 J/ R4 p) d\" }
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   8 S- o7 j\" a( j' W& i/ Q\" u
2. t=sys::clock();
   # M- P/ y. I2 s( L* s\" k4 e* N* L
3. s=0,x=0, & Y  Z% @6 d2 B' W& W5 P# Y
   
4. while{x<=1, //while循环算法；
   ; o# `9 A: C3 W* N' J
5.     y=1, 1 z  Q6 r5 o( N* T  w
   
6.     while{y<=2, 8 f6 Z5 }2 L0 x8 c\" g
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   * c: Z4 L- H1 r# w+ f# @! L
8.         y=y+0.0009
   ' N\" e! d$ V+ u- V! F1 f6 P
9.     }, # Z. y/ B, d1 h/ u5 `; C
   
10.     x=x+0.0011
    6 B7 }) B/ r9 D, q0 c3 e8 S+ f: W4 a
11. }, ) s) @! _\" D! B1 F7 \0 Y: `9 i
    
12. s;
      Q# t7 o/ a5 J' f
13. [sys::clock()-t]/1000;

复制代码

结果：
* ]) h' ~$ @2 h: I1008606.64947441
4 g/ ?% n0 P5 J! Z) w* G0.734   //时间，秒

0 d* h8 Y' U( G3 N我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

( F. _/ Q0 d2 k4 h' S/ J* B－－－－－－－
0 d- M2 }1 n0 N0 n# S5 ?( e& R& Z1 o
# G+ X+ [8 R6 \& z5 jForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   6 \: e6 R. G$ q# i' [
2. t=sys::clock();7 @# Q4 b1 i; t& Y
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   , u2 ?\" I' E% t( ~  l
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   + U1 a$ }0 K& V6 r9 @+ H
5. [sys::clock()-t]/1000;

复制代码

结果：
# Z: x; t8 |9 E$ ~% w  ^5 K1008606.64947441
3 F9 ~, K0 J1 o! x# N# ?* x  b: G0.719   //时间，秒

forcal

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

' R$ b6 q9 G1 G. |' B/ Q' G7 o2 \. Amatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! J0 q# O+ T: ]1 W3 y
   
2. tic
   ! W\" ^& S, d. n/ ?+ [- o. Q
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);& u& ^9 P8 }5 Z9 O3 S2 o
   
4. sum(sum(arrayfun(f,x,y)))
   ) ?' Y- D) u/ ~\" P( G' M
5. toc
   - Q% T1 k5 G) ]# K
6. & Q# e) D/ {) B( M$ x8 O
   
7. ans =
   5 `/ ]# K7 I, V0 j
8. 
   $ n* \) L# k+ v
9.   1.0086e+006
   & @0 x$ e% V) o5 s+ x3 N
10. + J5 R) {+ W, I; \- a
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];: Z: @/ x6 g3 s4 m8 |
   
4. mvar:! E1 O: W5 \) H
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. - u4 r: x2 U/ W% D
8. t=clock(),- [8 A6 I9 k: w) p2 j+ T) G5 A
   
10. oo{
    
11. % C5 W/ x0 `& K7 P4 _
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],7 b  U* H- |1 Y+ w
    
14.     arrayfun[HFor("f"),x,y].Sum[]3 l8 Y2 i6 ^# v0 T# [( C! Q4 d# i  l! z
    
16. };
    
17. / V! d+ C' Y6 p; v* p
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">: Z: @/ x6 g3 s4 m8 |</font><br /> mvar:<font class=\\"jammer\\">! E1 O: W5 \) H</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\\">- u4 r: x2 U/ W% D</span> t=clock(),<font class=\\"jammer\\">- [8 A6 I9 k: w) p2 j+ T) G5 A</font><br /> oo{<br /> <span style=\\"display:none\\">% C5 W/ x0 `& K7 P4 _</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">7 b  U* H- |1 Y+ w</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">3 l8 Y2 i6 ^# v0 T# [( C! Q4 d# i  l! z</font><br /> };<br /> <span style=\\"display:none\\">/ V! d+ C' Y6 p; v* p</span> [clock()-t]/1000;

结果：
: G( R+ e) V5 w1 ]2 b1008606.64947441
* c/ H' I0 f# l# V0 o0.735  秒
2 ]8 O7 a, L# [+ Z( |
! V; R; M  O6 t可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
. l' O8 @; R5 ^! T3 J
* e9 g' E) U! ?# h－－－－－－－－
+ J5 l* G9 W% I+ d  @$ a
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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