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

forcal

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

1. //用C++代码描述为：
   5 `* ~7 q- S2 s0 V
2. s=0.0; ; I3 K5 T  t8 X# I/ p- C
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   5 E$ {4 s+ K, ~& j9 Q
4. {
   3 }' V. e\" f+ e$ r
5.    for(y=1.0;y<=2.0;y=y+0.0009)# e( ]/ ?) {) ~0 c3 w( u% P
   
6.    {
   . W% o$ x; X$ d& H5 l/ u
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! m) L% g7 @& Q! r, l& j! j7 M
8.    }
   8 R9 n\" T, x& ~5 _
9. }

复制代码

Matlab代码：

1. tic
   ! ^+ h' ?: T- u6 O# @5 E2 @$ r: S. j6 S7 {
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);/ |* o7 A9 s$ s
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))& U8 n. u/ r4 F1 a9 c) \. H
   
4. toc
   3 u5 i0 d$ ^- p+ E: A
5. : n- o6 x5 g) X$ X
   
6. s =7 z3 z; c' z- w\" A: i# r
   
7. 0 ^5 C+ A* M# ^2 ?: D$ }. p, s4 V
   
8.   1.0086e+006* [! ]2 B  j0 C! J
   
9. 
   3 C8 B- X. [( Q7 L0 D1 {) o3 [
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];5 d\\" u( y+ p$ D4 M2 Y2 ?
   
4. mvar:7 @& b* q7 v& u0 u( e
   
6. t=clock(),
   
7. ) W5 D( r9 ?' l4 `0 a8 K+ Q
8. oo{+ r& |& h\\" N) r  {% L) E
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 7 |( y- d6 d0 c4 h
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]# d9 p! |1 j. r4 h1 }
    
14. };
    
15. ! B1 B+ a  N7 R  M0 {3 ]
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">5 d\\" u( y+ p$ D4 M2 Y2 ?</font><br /> mvar:<font class=\\"jammer\\">7 @& b* q7 v& u0 u( e</font><br /> t=clock(),<br /> <span style=\\"display:none\\">) W5 D( r9 ?' l4 `0 a8 K+ Q</span> oo{<font class=\\"jammer\\">+ r& |& h\\" N) r  {% L) E</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">7 |( y- d6 d0 c4 h</span>     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\\"># d9 p! |1 j. r4 h1 }</font><br /> };<br /> <span style=\\"display:none\\">! B1 B+ a  N7 R  M0 {3 ]</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641

' {" a3 E$ n5 ZForcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。

3 m, `6 n4 @, j* d5 v* aMatlab代码：

1. tic+ q3 h' E* m$ L$ F* ]% F
   
2. s=0;
   $ G3 v' Z9 t7 ~# D4 E; {
3. x = 0;
   + g- N  J; v/ e
4. y = 1;7 ?8 k\" `& e0 |! {/ r1 h  `
   
5. while x<1    ! z0 K0 Q; O4 M* q* z. a' D8 ^
   
6.     while y<2;        
   # ]! X! Q' }9 R7 j( l
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  t; e+ o3 s$ e
   
8.         y = y+0.0009;        3 H: F; I/ U) c1 x4 w! ~3 K, k
   
9.     end/ N5 j) y7 k) I& i
   
10.     x = x+0.0011;5 B8 [- X6 l1 ^0 B5 e1 h
    
11.     y = 1;( \8 s0 O7 h# U8 o; j6 V8 |, Y3 q- O
    
12. end5 k) l6 G4 B) Y. U9 J; ~
    
13. s1 m4 L- j! b; D3 }' v8 M
    
14. toc
    - w1 o2 v! L& I! [) p
15. ; C% W# Y' j6 U- f% c4 ?) ]8 d8 n4 b9 k
    
16. s =& s; I5 m' J- Y1 M9 M
    
17. ( V& L  ~( q: G( p: u
    
18.   1.0086e+006
    # ~# y) d& N# I: V9 q
19. , s' O\" O1 j2 E- x! J1 e
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   * u1 }) x8 c6 O$ z) u0 q
2. t=sys::clock();2 B9 U, Q\" l\" X# S
   
3. s=0,x=0, . e# V. l; N- F/ p1 i# x
   
4. while{x<=1, //while循环算法； ) x9 f# \: }  y\" C' q\" M$ K
   
5.     y=1, 6 }3 }4 G. t\" j5 O' a6 F
   
6.     while{y<=2,
   ; ^6 T* D6 K4 f% d+ t  D6 J
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), # @\" A# k! n5 B; n! _0 e
   
8.         y=y+0.0009
   4 N; N9 l; T8 [1 @5 V! q  z
9.     },
   + J' R% o/ Q, d1 h+ `
10.     x=x+0.0011 * W7 a* j/ j) s% i0 T# A
    
11. },
    4 G  l: K$ n' U. Y- |' }) w
12. s;, h9 N+ Y) v3 d* [  p
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
0 x: Q3 W( {0 {* @3 [
) @1 g3 W+ @- X4 k我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
& H9 f( \: E4 `) U. \6 ?4 Y$ T
( O* c( D5 N, x; m* P4 Q1 G  s. d1 O0 C－－－－－－－

5 Q2 b3 I: o& N( Q8 ^! J1 d2 ]5 AForcal中还有一个函数sum专门进行这种计算：

1. mvar:  H! l! d; z, x! `9 k3 V% y
   
2. t=sys::clock();
   / E1 t; Z) B( W( \
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); , H4 l\" K' k3 c) k! |) _( g6 A\" Q
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];+ |. V* z! Y; h' n  H* D2 g
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1 F2 f! ?* O6 R+ O! f8 p1008606.64947441
6 A' u4 v6 @2 ]% ~9 S9 j. G0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
" ?4 p! {& y6 R$ X
4 \( [6 p) k! ], S; Wmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));9 J9 Q6 Y+ X3 X4 h0 H
   
2. tic# p& O, v6 f& ]5 ^
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ' p' }( u* h1 F' D0 |
4. sum(sum(arrayfun(f,x,y)))
   / w2 {* V; S3 ?+ i
5. toc) }0 ?. O; f\" f2 C\" K7 V
   
6. . m- K0 M0 ^6 s% f1 P
   
7. ans =
   . w5 `1 _% `5 d  X
8. 
   + ?% S, ]/ q* a& f9 W
9.   1.0086e+0066 N7 f/ D: l  T, j2 U- R( f
   
10. 
    6 y0 A\" M' h1 G& W\" u
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3.   ]6 {8 i5 [  S
4. mvar:0 m\\" H& F4 `\\" S  H0 U, ~2 c
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); \\" u) H8 m+ f7 C0 J' Q5 i: u
   
8. t=clock(),
   
9. 4 g3 @( `$ R4 ~5 R$ v4 V3 N7 T
10. oo{
    
11. 3 G. p) F) b  ~4 [
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. % T( a  s% _\\" k4 S9 Z( e
14.     arrayfun[HFor("f"),x,y].Sum[]8 _5 j! C. t% R, I& f1 s8 }4 B
    
16. };+ X9 v) _\\" _4 k
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">  ]6 {8 i5 [  S</span> mvar:<font class=\\"jammer\\">0 m\\" H& F4 `\\" S  H0 U, ~2 c</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">\\" u) H8 m+ f7 C0 J' Q5 i: u</font><br /> t=clock(),<br /> <span style=\\"display:none\\">4 g3 @( `$ R4 ~5 R$ v4 V3 N7 T</span> oo{<br /> <span style=\\"display:none\\">3 G. p) F) b  ~4 [</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">% T( a  s% _\\" k4 S9 Z( e</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">8 _5 j! C. t% R, I& f1 s8 }4 B</font><br /> };<font class=\\"jammer\\">+ X9 v) _\\" _4 k</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
$ F# j# ~  e3 g) k5 _& g# N- D
$ n- y- m2 L8 i- O8 _5 I+ d可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
$ [$ t# x5 W  v; L4 Y8 a
  `8 v2 k# F9 l% I$ e－－－－－－－－
3 P' g& q& @9 j3 _
+ D* E3 r- Z$ a从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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