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

forcal

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

1. //用C++代码描述为：9 ?# n/ B, Y  I1 _$ M
   
2. s=0.0; ) ?) w  }* t4 U7 y8 v/ i
   
3. for(x=0.0;x<=1.0;x=x+0.0011) 6 V0 P8 e6 `) {4 P
   
4. {
   * W2 o+ ~1 K3 z
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   % A1 p% Z- _8 U4 f7 i, H
6.    {
   7 P3 o8 h/ ~2 E! D
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));' A- m1 Q: U4 B3 N8 `6 U+ T: |( o0 @
   
8.    }  R: j# ~\" W; N0 A6 R: ^
   
9. }

复制代码

Matlab代码：

1. tic
   8 l3 P  s4 }- |6 \9 }
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   7 T7 L$ Z, {9 D8 j& F
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))). Y. Z+ `9 s  `
   
4. toc7 M2 K+ H& }$ u4 [9 j9 w& _
   
5. ' h6 a& i% H1 H8 ]5 E3 G
   
6. s =5 l. @. Z1 Y) l7 g/ k
   
7. 
   9 |1 k9 E/ ]! e: ^
8.   1.0086e+006+ ^4 b* M0 {. ^, ]! c: A
   
9. 
   5 [\" x  ?: m, }; F8 Z; q8 B
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. + s0 N( r- z  m
4. mvar:7 ]0 @9 O. i; f$ f
   
6. t=clock(),
   
7. + W3 z0 Q) _% b
8. oo{
   
9. + o% m2 j& E# \9 V
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],4 ~; a# V! r( v* U& y$ j# g
    
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. ; ]4 N+ S  z# m2 R
14. };
    
15. 0 w/ c# F4 z- q: {( ^\\" x
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">+ s0 N( r- z  m</span> mvar:<font class=\\"jammer\\">7 ]0 @9 O. i; f$ f</font><br /> t=clock(),<br /> <span style=\\"display:none\\">+ W3 z0 Q) _% b</span> oo{<br /> <span style=\\"display:none\\">+ o% m2 j& E# \9 V</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">4 ~; a# V! r( v* U& y$ j# g</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\\">; ]4 N+ S  z# m2 R</span> };<br /> <span style=\\"display:none\\">0 w/ c# F4 z- q: {( ^\\" x</span> [clock()-t]/1000;

结果：
1008606.64947441
5 n" r/ _7 A( h5 C+ c' }0.641
8 B1 L& f- s( q2 H8 _; |+ K
Forcal比Matlab稍慢些。
% C( ]4 c! }8 ]* _
5 |/ }% U$ }0 i! w" D2 Q! E2 x－－－－－－－－－－
+ V$ @" a4 w' G
再看循环效率。

Matlab代码：

1. tic5 E: ?1 _2 K% }- ~
   
2. s=0;
   - B2 d\" M+ j# T% M
3. x = 0;/ |. J1 U( P; G/ Y' y* n; j
   
4. y = 1;! P; a1 X5 {1 g6 p+ t
   
5. while x<1    % o  x4 K% S+ I6 i3 O! m  B
   
6.     while y<2;        ; _9 E/ H9 H\" a4 {
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));; P* a9 p7 E3 q& J6 W& s4 x
   
8.         y = y+0.0009;        
   , h& m8 u3 M& m* ~; F+ |1 R0 V
9.     end
   0 N% ^) K: f4 @% `
10.     x = x+0.0011;
    5 |  J' t0 ]# Y
11.     y = 1;' H* j3 X8 _# I, c0 s
    
12. end% `. d8 T7 B' |! \, s4 L, u
    
13. s
    5 M, b6 u; N! I( b+ V9 J
14. toc8 V' A2 L8 }6 ~8 D, ^5 s
    
15. 
    0 V  U4 ?+ Q- m
16. s =
    : |2 k. X! O- ]; y% g4 |5 ~
17. # x- @$ n  \! Z: T4 p; ?; ~
    
18.   1.0086e+006
    6 K7 j' K$ x. T# g# b4 T# s\" i) q2 J3 K
19. : k) g: W& ?& ]/ e1 ^4 L! `3 }
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:- C$ [5 R5 N& v0 C8 H\" u) c
   
2. t=sys::clock();' N) H# |/ Q* u- g# J
   
3. s=0,x=0, 4 G/ N; d; c  C7 s, ~( `5 ?/ W) S
   
4. while{x<=1, //while循环算法； * C5 G( F3 O* h. M+ I
   
5.     y=1, ' A( G/ L5 {4 N3 h1 z. f' S
   
6.     while{y<=2, - I5 l- [  a5 s# y
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 6 B; p$ ?8 D: |- d$ G% j$ r
   
8.         y=y+0.0009
   . M0 g: E) K2 I* V- Y1 g* I5 R  k
9.     },
   8 f' d/ ]1 W/ u
10.     x=x+0.0011
    5 ~1 m/ K3 `4 M; T- C
11. },
    ; y5 K\" k3 [0 u# ~+ \
12. s;3 t8 @7 t* w6 e, W8 w6 n2 k( I9 _
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
- B3 n: f5 ^3 `- U. h4 }3 x* S
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

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

1. mvar:- n6 l' {$ Y  W5 D
   
2. t=sys::clock();
   - t% C% k1 ~; J4 B& d4 d' z7 }: R3 Q4 R
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
     X7 I( E- e* C& P# I/ [  Y
4. sum["f",0,1,0.0011 : 1,2,0.0009];\" s# l& |. a- C7 B- B
   
5. [sys::clock()-t]/1000;

复制代码

结果：
) f# m5 F2 v1 T1008606.64947441
0.719   //时间，秒

forcal

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

6 I+ [' f  E( q% n$ Z& hmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));\" P( ?: _( Q# E) e7 B6 {
   
2. tic; a6 p3 Y1 v/ z6 s
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 U/ [5 e2 |; W9 F' O& J5 k( t. i
   
4. sum(sum(arrayfun(f,x,y)))
   # Z: M+ y6 j8 p5 z! y5 J7 O* R
5. toc
   & n8 v0 W4 m' p' ]1 s7 y
6. 5 v% a. W! _% x2 U
   
7. ans =9 C\" Y9 z3 y$ ~$ W
   
8. 
   ) V3 q( u/ d, L% D5 W+ `
9.   1.0086e+006
   ! N* g6 B* s. n/ h2 \! q
10. 
    % F) s5 M6 c% z/ c& \  [& V
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 2 {( M$ V4 x. n% {+ i; K2 {' |
4. mvar:
   
5. % p% W/ {+ k  Z( f+ q: N5 x: k
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. ; h2 a8 w6 q5 X/ L1 t
8. t=clock(),
   
9. - U! h$ p5 [- Z0 r9 Q, f
10. oo{
    
11. : Q# ^6 j+ N0 u# h. x
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],8 s0 _% y3 @/ I, S/ K! B
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. - s! c7 b\\" u) a, C
16. };' K3 `' }; J8 B3 c
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">2 {( M$ V4 x. n% {+ i; K2 {' |</span> mvar:<br /> <span style=\\"display:none\\">% p% W/ {+ k  Z( f+ q: N5 x: k</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">; h2 a8 w6 q5 X/ L1 t</span> t=clock(),<br /> <span style=\\"display:none\\">- U! h$ p5 [- Z0 r9 Q, f</span> oo{<br /> <span style=\\"display:none\\">: Q# ^6 j+ N0 u# h. x</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">8 s0 _% y3 @/ I, S/ K! B</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">- s! c7 b\\" u) a, C</span> };<font class=\\"jammer\\">' K3 `' }; J8 B3 c</font><br /> [clock()-t]/1000;

结果：
! X7 d4 t. v0 B, j0 f! h1008606.64947441
0.735  秒
! x2 l# T3 |- t: z
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
9 E/ z0 e. _$ n0 I% t0 N( o7 c' D
－－－－－－－－

- n( G2 w: p3 x& g2 r从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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