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

forcal

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

1. //用C++代码描述为：' p% B4 t* ?- p; |- w2 h\" O
   
2. s=0.0;
   + }\" E) D8 j& t6 G
3. for(x=0.0;x<=1.0;x=x+0.0011) ( B: M\" _' ~& Q2 T! X$ y0 r
   
4. {
   2 [. L2 N$ \% Y8 l- D& K$ H
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   ! {  m# F; S4 a
6.    {  l1 Q9 A) A$ H% v2 w/ M0 g
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ n2 `' D9 e7 J2 i! k, j5 k
   
8.    }
   ( e( B) Q: w: }3 T1 I6 T8 j
9. }

复制代码

Matlab代码：

1. tic
   $ D9 Q8 Q2 e% q
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   \" ?1 f$ j6 o$ P* x' K
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
     C* i7 K/ z8 |/ Z. e8 L3 C  i. g
4. toc+ h' r\" O1 ]$ x- v7 {
   
5.   B3 Q; B, F, i( N7 O* j7 N7 V
   
6. s =* T% I7 g0 t( z2 ^) A: ]- D/ @9 O
   
7. 
   ( E3 U' W5 t\" m
8.   1.0086e+006/ E  Z* j4 s) N0 L\" m/ N5 i1 s' \\" b
   
9. 
   / y  L  ^/ W  `, n$ d, L( |# c
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 8 e4 _# w\\" g4 ]# w
4. mvar:6 B( K. w4 d: {1 }4 r
   
6. t=clock(),
   
7. + z1 L, [5 I2 p  i
8. oo{
   
9. 5 B% O) l0 E* _2 r; s
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],% v/ N! W$ T- Y! Y
    
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]' Q3 _$ U% r4 N1 H% w
    
14. };: f\\" v4 c1 X0 `5 P8 L& Z; T3 e5 ~
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">8 e4 _# w\\" g4 ]# w</span> mvar:<font class=\\"jammer\\">6 B( K. w4 d: {1 }4 r</font><br /> t=clock(),<br /> <span style=\\"display:none\\">+ z1 L, [5 I2 p  i</span> oo{<br /> <span style=\\"display:none\\">5 B% O) l0 E* _2 r; s</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">% v/ N! W$ T- Y! Y</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\\">' Q3 _$ U% r4 N1 H% w</font><br /> };<font class=\\"jammer\\">: f\\" v4 c1 X0 `5 P8 L& Z; T3 e5 ~</font><br /> [clock()-t]/1000;

结果：
- e; l& Y( }& i" ^) A8 I0 D1008606.64947441
7 n" V' n0 G9 g- G* _8 v; D' P0.641
. T+ T3 w# O+ R% J3 D- U
% A- X% z/ d4 U0 N2 T: F& pForcal比Matlab稍慢些。

－－－－－－－－－－
5 v# n  L' X1 t% N1 n1 o2 \: }+ H
再看循环效率。

- S+ w7 y0 d8 f! _# V+ MMatlab代码：

1. tic
   ! j3 |6 h, j1 j/ i& \
2. s=0;3 X0 V! O+ [; \: d* j
   
3. x = 0;' }\" S  \' Z' s+ H$ u' s% T: L
   
4. y = 1;4 v$ D6 L, o) c/ t3 C/ q\" @% k
   
5. while x<1   
   0 f% S3 @+ h) X% F\" U$ ?1 H' h
6.     while y<2;        
   / e( n7 b% W' ]/ M! ~+ Z
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ e, b; m% T8 Z$ I8 B
   
8.         y = y+0.0009;        
   & m+ V. b\" f) Z7 }# N  X, `# R0 w
9.     end* F! z: Z. m8 V1 y
   
10.     x = x+0.0011;& o) g/ d! P% b\" w7 a5 k7 ^
    
11.     y = 1;
    \" \- O' D6 X\" N% ~1 C3 Z
12. end) k% {9 y' @; I; c4 \8 P
    
13. s6 S& [1 r/ Q+ O  @# z* E# @
    
14. toc+ ]& d\" h* ]( @3 ^4 z
    
15. % J\" a- ?/ K5 M4 e9 J3 E
    
16. s =
    ! i7 y' r- B/ I
17. 
    4 s  D8 q0 u: S\" y' \
18.   1.0086e+006
    . w; F. g# R\" k, X
19. 4 c5 D. B2 t\" r$ Z
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:# i1 J. O+ L( H3 E) ?& ?# ]5 a3 n
   
2. t=sys::clock();
   & D5 J; g+ T* i1 N
3. s=0,x=0,
   ( E5 u6 ]8 h5 A. C
4. while{x<=1, //while循环算法；
   0 Z# m# h& V# R! j
5.     y=1, + v0 A3 M& O& S% i
   
6.     while{y<=2, / R6 A, J) L. u, G6 O
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), . L, ^$ c- s/ q- `( w- M0 O5 W
   
8.         y=y+0.0009 + g/ l0 A# w( B: }% K- Q9 q
   
9.     },
   4 }9 p( ?! e. Z& D
10.     x=x+0.0011
    ! Q3 D5 D& c( P
11. },
    ( d2 e6 }! Q/ \, h1 z/ B
12. s;
    ' u5 c+ ~' Y5 h+ o- X& V
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
& m& f2 N/ [' @9 ]
! p( t: i! ?( A) w' [$ j, B6 N我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

* h- n) t' D2 f－－－－－－－

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

1. mvar:
   + G, u, ^\" K4 d
2. t=sys::clock();
   5 t& S, a7 w! I# O
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / C4 z# d4 s8 e
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   : A' f$ Z6 B; ]% Z8 K
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

2 z: }; J  S- }, w- nmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % J- |' ?' }\" \5 D& x4 k0 Y
2. tic# k9 ]* k  @$ V/ Q9 b7 I- n
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);  p5 S: D! Y- Y9 M8 `, ?
   
4. sum(sum(arrayfun(f,x,y)))9 `6 ?: V9 R7 C( q
   
5. toc
   , J* D3 z$ J* l- {; z\" J1 i
6. 9 l. }$ c6 ^  P. m
   
7. ans =
   - W7 Q. v6 A5 l# K( c6 d6 [
8. 
   ! U) Z5 B. V; w, f
9.   1.0086e+006/ u4 ?$ A5 v) A, c. l
   
10. 
      T) V7 l+ X0 x0 F\" [; Z9 O  J
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. , t- P% o0 J# \
4. mvar:
   
5. - O1 |! ^9 L& s
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); \\" j  ]' G: y1 Z$ R
   
8. t=clock(),\\" O7 V7 L$ B& ~6 O1 b& o  H: h9 H
   
10. oo{
    
11. ! z! p& k! x% ~& O& ^1 i2 t
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. & n\\" w\\" ]- r0 g4 ^) h, F3 G$ T
14.     arrayfun[HFor("f"),x,y].Sum[]9 Y\\" p7 J) ?4 }; u' M! M
    
16. };, W3 |6 s: p' u
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">, t- P% o0 J# \</span> mvar:<br /> <span style=\\"display:none\\">- O1 |! ^9 L& s</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">\\" j  ]' G: y1 Z$ R</font><br /> t=clock(),<font class=\\"jammer\\">\\" O7 V7 L$ B& ~6 O1 b& o  H: h9 H</font><br /> oo{<br /> <span style=\\"display:none\\">! z! p& k! x% ~& O& ^1 i2 t</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">& n\\" w\\" ]- r0 g4 ^) h, F3 G$ T</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">9 Y\\" p7 J) ?4 }; u' M! M</font><br /> };<font class=\\"jammer\\">, W3 |6 s: p' u</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0 Z# w: Q& _+ K" ~0.735  秒

6 ]; P, H" Z: L1 T" c! [( L3 E! l可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
. y7 z1 X  d# b( M0 Z2 N
－－－－－－－－

从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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