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

forcal

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

1. //用C++代码描述为：% A+ c4 g  ^  w
   
2. s=0.0; * u* n) q2 E7 X& l% z
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   4 T, I2 U0 P1 c7 l2 i* z
4. {! k$ h7 C6 \6 c( G0 k5 [. ]
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)3 V$ J) V* W! r5 _  e6 J
   
6.    {
     ?$ P' V: `9 i% a* ?
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 @1 J, d0 y2 K+ R7 _1 g' m9 F
8.    }\" h  F# x- u( ]4 @2 j6 W$ j
   
9. }

复制代码

Matlab代码：

1. tic) Y' ~\" G* m% B- @. V
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   1 x+ h4 T4 L* y: g) b8 b& F
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))! `7 H% m1 y% f! q+ @7 q4 d
   
4. toc
   . G\" y* U  _! c1 Q
5. 9 i' e% G8 e1 q' D
   
6. s =
   ( C* W6 Q8 b4 Y+ L  `3 n
7. 
     n* F) x7 Y& N( i5 ^- b* L, ~
8.   1.0086e+006
   * m8 c$ c2 o8 x$ E$ o8 r
9. 
   . q8 y/ z! }9 k! ]# v
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];2 D5 e% o3 }$ o. Q( T# S8 N
   
4. mvar:6 ?4 j6 G4 P. V$ l% C
   
6. t=clock(),
   
7. - E7 g\\" m8 a! H. k) u% }
8. oo{
   
9. * R& T- ?3 ?+ W
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],! f5 ~+ x3 b3 A. X7 Z7 I
    
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. ; w9 q' R6 H  V: ?
14. };! y8 k% V6 v$ ]7 a% _: Y
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">2 D5 e% o3 }$ o. Q( T# S8 N</font><br /> mvar:<font class=\\"jammer\\">6 ?4 j6 G4 P. V$ l% C</font><br /> t=clock(),<br /> <span style=\\"display:none\\">- E7 g\\" m8 a! H. k) u% }</span> oo{<br /> <span style=\\"display:none\\">* R& T- ?3 ?+ W</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">! f5 ~+ x3 b3 A. X7 Z7 I</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\\">; w9 q' R6 H  V: ?</span> };<font class=\\"jammer\\">! y8 k% V6 v$ ]7 a% _: Y</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
5 i. Y3 z! f+ J1 O: Q8 z9 U0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－
& q- ]. H7 Y5 Q: u) r1 H5 d
再看循环效率。
. t- w  G! o" j$ p1 E5 t& I
# O. x3 i$ d. g! g1 s/ LMatlab代码：

1. tic
     Y1 o1 j6 x  [1 L1 [) ]\" Q
2. s=0;
   5 w+ U- _/ k8 M: J
3. x = 0;4 A5 E4 e% I9 L4 T9 q2 ]
   
4. y = 1;, Z3 O6 r/ f$ e5 H; ~
   
5. while x<1   
   % x3 c1 I& T, S/ K
6.     while y<2;        4 n+ F( x0 d\" }
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   # f\" t) L  t, u& |1 _# z( A: R
8.         y = y+0.0009;        
   ( p3 o3 f\" c9 H$ W7 q- a+ G
9.     end
   4 s1 v, M+ H- t0 ]
10.     x = x+0.0011;6 M% ?+ M% A$ c% K6 W5 N+ d) r
    
11.     y = 1;# G+ a2 R+ m: ^, W. `' y
    
12. end
    . W9 E- C- n  ~- x
13. s+ _$ T. k& @\" S\" E' q
    
14. toc
    $ ~% V' b' m! O0 l- G+ B
15. ) s( J$ G* Q0 r! [% t
    
16. s =! ^2 G- P, u1 Z' k
    
17. ; Y+ v2 j8 v  h: {
    
18.   1.0086e+006
    9 v/ Q) _/ y6 R\" v. u
19. 
    $ G' A1 k- Y  Q; o3 O2 O/ F
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:: f1 k& v) F6 L; Y* e
   
2. t=sys::clock();
   \" p6 M1 E2 p  W
3. s=0,x=0,
   ) w' X' |4 V5 \) h# n0 X: ?
4. while{x<=1, //while循环算法；
   ( c# H- R4 [\" ?4 j3 t9 v. K
5.     y=1,
   , M. S& i* j+ h6 m+ n( N
6.     while{y<=2, 4 ]( p( m0 j3 A3 M\" f  V
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ! l) I5 H* I& P. z4 X% u
8.         y=y+0.0009 % x3 n3 x; C\" O2 ~' S5 X
   
9.     }, % v! U8 ]' o( }! P! R/ z) I
   
10.     x=x+0.0011 / S9 C# Z1 v4 `
    
11. }, # U- q# l! {; @, j6 s
    
12. s;
    ! ~' ^: U1 i. W6 {; [9 f/ L
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
: B8 ~6 @* R+ z7 T' H; S0.734   //时间，秒

/ l: m5 m# f# y, ~; N: s我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
* t$ H( J# V8 k; f* N$ ~0 O
+ o: \7 l- L* J3 f- v$ D9 gForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   ) }- s! v* D( a+ {' E\" {
2. t=sys::clock();
   7 U% E! k5 D! }/ n
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . p# f: ^# {0 Z9 j
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ; V: k+ }$ `! M. X# Z
5. [sys::clock()-t]/1000;

复制代码

结果：
" p: E6 C3 |0 V9 T9 _1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
4 n2 D$ V2 q/ F" o
; n, z1 `3 h- ]8 f  Gmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / w( W: d$ Q5 H
2. tic# }! D0 L; c9 w; L4 V, T
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   5 P' g2 @2 E! w$ @
4. sum(sum(arrayfun(f,x,y))); a8 e( N' h  P\" x* M) C
   
5. toc$ }( ?7 t1 F( H
   
6. 
   6 c$ p( k; ^2 ^. r. g
7. ans =  C. }2 X' t  E+ P* @
   
8. 
   ; t% l, `5 n! k% E/ ~
9.   1.0086e+006
   5 S# ]1 N# J8 k( l+ I! G
10. \" l; `2 ^  K- j
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. % ~. B\\" R) L  c3 M
4. mvar:8 r. M/ r9 O$ w7 V2 v: N
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 2 ^- D; T1 g. d7 ?& z
   
8. t=clock(),
   
9. $ ?* i& N0 E+ |1 Z
10. oo{
    
11. ( R/ U% y4 F- F( I4 B, J
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 8 M1 c\\" D! o& r) L$ p: b
14.     arrayfun[HFor("f"),x,y].Sum[]7 K\\" I$ F, ~3 y5 e9 j; S
    
16. };
    
17. - v( K: W6 E3 C& `; ~
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">% ~. B\\" R) L  c3 M</span> mvar:<font class=\\"jammer\\">8 r. M/ r9 O$ w7 V2 v: N</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">2 ^- D; T1 g. d7 ?& z</font><br /> t=clock(),<br /> <span style=\\"display:none\\">$ ?* i& N0 E+ |1 Z</span> oo{<br /> <span style=\\"display:none\\">( R/ U% y4 F- F( I4 B, J</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">8 M1 c\\" D! o& r) L$ p: b</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">7 K\\" I$ F, ~3 y5 e9 j; S</font><br /> };<br /> <span style=\\"display:none\\">- v( K: W6 E3 C& `; ~</span> [clock()-t]/1000;

结果：
1008606.64947441
+ d( Q8 e, |& b& W0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
2 l0 p$ {$ [) V: `/ _
－－－－－－－－

4 H, p0 s2 l- w# B) r0 e从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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