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

forcal

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

1. //用C++代码描述为：
   6 V  v& U. d& D
2. s=0.0;
   & O' @% W# Y* O! P5 w
3. for(x=0.0;x<=1.0;x=x+0.0011) 2 F6 \, N5 ^. m- a* R
   
4. {
   \" V/ ]7 d6 V4 o# I! o1 b
5.    for(y=1.0;y<=2.0;y=y+0.0009)) h& }\" J8 ^' w1 U! y
   
6.    {
   % y: z$ A1 a9 U! Z& r; S1 Q\" G1 b. U
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! h# F5 \) h3 C0 D$ X- ]
   
8.    }
   + k& j, b) e( o' z2 ]% k
9. }

复制代码

Matlab代码：

1. tic
   : L6 Q8 `4 {: S! q; m& c
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * H# }' e8 t9 O1 \- d& _
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   8 }+ r/ n2 n$ S* U
4. toc
   / C% ~. j* Y, n. |5 l
5. : y# B# Y# Y$ L
   
6. s =+ e0 ?9 K4 U1 x- _- _2 t
   
7. ; T, k1 c8 J* K+ I/ P
   
8.   1.0086e+006
   ; r\" Y! E( U# F
9. 
   3 H+ L. Z/ Y! J
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ' e\\" ^\\" K$ ^0 }6 C
4. mvar:+ z; ~! I- M& C7 s! o5 r8 h4 I
   
6. t=clock(),2 F! {' b' P$ @% @
   
8. oo{
   
9. ; n* K/ z* a5 a9 j4 I$ ]: V
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],- C( t/ X6 m: l, P# M# e
    
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. 9 l, \3 B+ U% [( t/ m* R
14. };& v, W+ E\\" u+ t( {
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">' e\\" ^\\" K$ ^0 }6 C</span> mvar:<font class=\\"jammer\\">+ z; ~! I- M& C7 s! o5 r8 h4 I</font><br /> t=clock(),<font class=\\"jammer\\">2 F! {' b' P$ @% @</font><br /> oo{<br /> <span style=\\"display:none\\">; n* K/ z* a5 a9 j4 I$ ]: V</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">- C( t/ X6 m: l, P# M# e</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\\">9 l, \3 B+ U% [( t/ m* R</span> };<font class=\\"jammer\\">& v, W+ E\\" u+ t( {</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.641
8 h! _& |7 W* D/ B1 T
Forcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic
   9 p' I' m; c4 e5 S2 Q0 ^
2. s=0;& W7 R( H\" X2 C2 \7 H7 T! H
   
3. x = 0;: O- v7 F, A9 [, K9 E: F
   
4. y = 1;) m( j9 A4 r% J# u
   
5. while x<1    : T( w8 N  e' ?8 c! H
   
6.     while y<2;        4 F% d: g& a- A. O( ]
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));/ a( \5 K* B7 k' ?! v
   
8.         y = y+0.0009;        
   + w/ ^2 p2 H  i+ j5 _
9.     end% E2 l2 Y9 j\" H! d* d; N. z! U
   
10.     x = x+0.0011;: _2 t2 _8 y/ c. }9 e
    
11.     y = 1;7 }) B, a, W5 Y5 ]. X
    
12. end5 c+ D+ c5 G2 S& B$ z; a
    
13. s) b2 Y# |+ G( P
    
14. toc
    8 [0 P) H+ s# i
15. 3 K( b$ U- \$ F& Q, ?! k4 H. l% x
    
16. s =
    2 k  \* E! C* i, {7 c1 F
17. 
    2 A! V+ g8 Q: j) D9 X
18.   1.0086e+006
    % k' n8 d\" |/ g- b) D
19. & G. c  b+ d( I8 ~( `5 |; O3 N& I' H9 b# ]
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:& {1 O- q: C! D# W# ]1 |* c
   
2. t=sys::clock();
   - a0 y! n8 K% e7 d
3. s=0,x=0, % ]\" `7 {. R; |4 A! v4 a
   
4. while{x<=1, //while循环算法； 7 d6 q* M9 j5 E1 b( B\" V
   
5.     y=1, ) G& C7 s1 i, [, x9 P5 u& ~7 s' ]
   
6.     while{y<=2, ) N/ _- r: F+ q* D
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   % g' U  y7 U2 r, X4 _
8.         y=y+0.0009
   2 [0 {/ h+ b) \/ W
9.     },
   & x$ g  y0 `/ s7 r/ Z( h
10.     x=x+0.0011
    9 L! Y4 ~5 H: x* L# F1 y
11. },
    + ]9 R7 p+ W7 V! f
12. s;
    * U3 B, U) g7 F( ?3 J' t0 s
13. [sys::clock()-t]/1000;

复制代码

结果：
3 P: w) O! P6 m: {1008606.64947441
0.734   //时间，秒
2 {& H6 e5 U, |& e) h8 I
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
$ B8 _4 t( k9 p; P! i3 u
0 E# u$ P8 ~! Y: S1 d－－－－－－－
6 k0 v- W/ L+ |* x7 b' `
) t0 V3 |5 i6 B. k  BForcal中还有一个函数sum专门进行这种计算：

1. mvar:6 n% c. N% |- t1 }! @
   
2. t=sys::clock();
   , b  q% M\" q! u( M9 M% O
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   $ h; M( i5 K5 I
4. sum["f",0,1,0.0011 : 1,2,0.0009];) D* S+ N. a, A\" v7 `0 H
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
2 D4 H' |8 H6 @: y6 X0.719   //时间，秒

forcal

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

- F( g: V9 M7 X4 }  Dmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - q# v& h8 `& ^, [( D
2. tic/ p  _* ?' H# F1 \1 i- i
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   % X* o/ N) W- Q* b6 Q
4. sum(sum(arrayfun(f,x,y)))
   \" A2 ~9 n3 b6 U\" h) R$ R
5. toc
   0 x7 n2 h  ~5 J& S
6. / M  E+ w4 g% B5 T. G0 f& _
   
7. ans =
   2 [+ V( _- g% z
8. 
   . }: B% i  d5 K' A0 o\" Q
9.   1.0086e+006\" g$ \+ P4 p; m1 N
   
10. 
    + A4 c' P$ r) b2 `6 ^
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];* f; }, g7 F# j4 w! l
   
4. mvar:
   
5. 4 n( _- L3 O; ?3 x
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); . y+ T5 w2 z: U
   
8. t=clock(),: t# j6 q: o, f5 v+ n
   
10. oo{\\" e1 F3 r6 {, }! e* m# X2 R7 g
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],3 s1 I# ?2 {* T+ v* n1 o' F' W
    
14.     arrayfun[HFor("f"),x,y].Sum[]7 q. S: }( y7 V. ?$ I7 Q
    
16. };
    
17. * J1 K+ j0 i\\" ?, B# z9 Q1 l
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">* f; }, g7 F# j4 w! l</font><br /> mvar:<br /> <span style=\\"display:none\\">4 n( _- L3 O; ?3 x</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">. y+ T5 w2 z: U</font><br /> t=clock(),<font class=\\"jammer\\">: t# j6 q: o, f5 v+ n</font><br /> oo{<font class=\\"jammer\\">\\" e1 F3 r6 {, }! e* m# X2 R7 g</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">3 s1 I# ?2 {* T+ v* n1 o' F' W</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">7 q. S: }( y7 V. ?$ I7 Q</font><br /> };<br /> <span style=\\"display:none\\">* J1 K+ j0 i\\" ?, B# z9 Q1 l</span> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
( I0 \3 r' t: j( c" m
: F0 q$ v" y& ]- P. V* P可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
2 o& v# a1 m- l2 t( z
$ f/ h- _( l! S% ~# o从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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