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

forcal

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

1. //用C++代码描述为：5 Z  d; e; h' ^, S% _1 k
   
2. s=0.0; , s4 L. `% B, C\" u# k0 r0 M
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   - S# S% `  g% o, l1 Q9 a- ^
4. {% ]) q- c) B) A, [$ ]; w
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)4 |; J* F( X1 J# Q9 |* Y) v
   
6.    {: B' V0 ^9 W8 M$ p( \
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));- ~5 E# c) ~) `6 p% C( @
   
8.    }
   # ?* }; I: }& d$ N
9. }

复制代码

Matlab代码：

1. tic
   + n1 N7 g8 Q* E; \8 o5 |4 X
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);9 u$ w/ b/ @3 X- B
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))7 l& x6 E! M: F) v# r5 d
   
4. toc3 ?4 }% N\" Z: k# s  q( @
   
5. 
   2 h) M9 ^% g- N* X! D
6. s =% j. l7 w! r, R7 p
   
7. 
   ; D\" s* T0 s* S
8.   1.0086e+006; D# L; T$ C. X8 P6 G$ O\" W' C
   
9. - _' O& W& S6 ]9 C1 o& ~& \  e
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];! i0 {\\" E0 h  h0 x
   
4. mvar:
   
5. $ N3 _8 u% K1 V) u6 @5 S
6. t=clock(),
   
7. % |7 b- f2 K) Y. ?
8. oo{
   
9. ) I. X& y% R6 u. P# h3 I, s: x
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. / `+ i' \9 k9 x: T
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. \\" s. Z, k- Z5 L\\" c; |& ^
14. };
    
15. 9 V4 \; }9 ]0 Z- m, {5 d
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">! i0 {\\" E0 h  h0 x</font><br /> mvar:<br /> <span style=\\"display:none\\">$ N3 _8 u% K1 V) u6 @5 S</span> t=clock(),<br /> <span style=\\"display:none\\">% |7 b- f2 K) Y. ?</span> oo{<br /> <span style=\\"display:none\\">) I. X& y% R6 u. P# h3 I, s: x</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">/ `+ i' \9 k9 x: T</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]<br /> <span style=\\"display:none\\">\\" s. Z, k- Z5 L\\" c; |& ^</span> };<br /> <span style=\\"display:none\\">9 V4 \; }9 ]0 Z- m, {5 d</span> [clock()-t]/1000;

结果：
1008606.64947441
# f! y5 n: k! U! R( j; w; q/ L3 k0.641

- |7 T. z0 L$ z; j, i, ^& BForcal比Matlab稍慢些。
( t$ _/ g% n* R  `
8 x% i' Q* I' u2 Y* a－－－－－－－－－－
7 Q; R1 M$ c( i9 I: M1 x6 F. X
再看循环效率。
* a7 Q8 d- |4 r+ _* K
Matlab代码：

1. tic0 y& T/ ^. N' Y+ F9 q/ q; ~/ R
   
2. s=0;1 j0 ~& l# j* U% S$ @& v8 t/ p
   
3. x = 0;
   + m: b% x1 J; L/ p+ I$ W
4. y = 1;
   ! z: q( N% }9 W  U: o
5. while x<1    ' z  q0 k, p: ^( b/ M
   
6.     while y<2;        : ~1 y  l5 P3 r# \' [0 e
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));( ^, R8 w8 B( G/ t, Y# t
   
8.         y = y+0.0009;        ! G: {3 @5 `. Q- r+ C( @
   
9.     end
   6 w: U8 y5 {' w+ J/ s9 ?
10.     x = x+0.0011;
    1 {; N; y3 S1 v* @* r* E* E8 H; |; `
11.     y = 1;
    ! x  ]. D7 i* [- x\" k# Y7 S# [
12. end
    ) x; T$ A8 f! {3 D5 }1 j
13. s
    + F9 b3 P  q  A2 a8 E& C8 t
14. toc
    ' I4 [8 M\" i, H& A8 R
15. 0 C* q& O# v- z
    
16. s =6 p5 y, ^& p1 T5 u7 u- p. I5 U
    
17. 
    0 R  j# z; [; ]9 W; N% Q' B
18.   1.0086e+006/ a8 [9 H4 I) b7 \, ^
    
19. 
    ' j4 B0 S! x' n- h! s, G9 ?
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:8 {5 }. I1 ]3 e% c
   
2. t=sys::clock();1 o6 |4 _: @, q% F# x6 u8 Z8 J# j
   
3. s=0,x=0, / x) A6 ]7 |2 T! [6 [
   
4. while{x<=1, //while循环算法； ) |3 p9 b3 ^; G7 E
   
5.     y=1, % }( W. y% ^7 M
   
6.     while{y<=2, . x  N7 l2 M0 V! ?6 |$ k
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ' Y/ S5 D. n8 I& n5 y
   
8.         y=y+0.0009
   % i9 {; G. K% j1 ]
9.     }, 9 ?- d/ L( Y7 f/ c
   
10.     x=x+0.0011
      a6 a3 ?2 c! k5 J) T% w
11. }, 9 ?4 i7 a) t( O
    
12. s;. y. x5 Y+ k2 {! |0 a
    
13. [sys::clock()-t]/1000;

复制代码

结果：
: x5 j! T0 C* P& z3 R! K2 ]: z1008606.64947441
! l. j6 |# l) F0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
  S0 c( U$ O. G% H+ C! K9 L1 i
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:2 O. b! f- u$ t6 I( a! x  ~
   
2. t=sys::clock();
   ( ^1 ~8 r& g9 W, E+ S( R
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 F\" H7 M5 o; z
4. sum["f",0,1,0.0011 : 1,2,0.0009];
     y\" w$ |6 s! w8 \0 x
5. [sys::clock()-t]/1000;

复制代码

结果：
: I) ~. v6 N! U. |1 F1008606.64947441
0.719   //时间，秒

forcal

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

7 ~$ F0 D! ?( o7 [" _0 z: rmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));6 u4 E& Y8 n( f- A, B% M0 n2 R
   
2. tic! ]) Q! T$ G0 n. [' e! g
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ! f1 K+ ~4 V5 f! U; }
4. sum(sum(arrayfun(f,x,y)))
   ; P, u4 @& F$ z6 M) J* d
5. toc
   , q* M  {* l# P. a( i' M
6. 6 G5 z) S\" J\" Y4 ~' T+ t7 f  s. s/ P* v
   
7. ans =9 }& X; e8 D# L
   
8. 
   7 b7 ~' g( L' Z, n2 B) ?  l1 C
9.   1.0086e+006
   ( j, `+ p\" A* K3 s* K
10. # x+ `. U8 ?4 ~$ L% Y\" j
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];4 J! s: Q1 ^& M6 h
   
4. mvar:% [4 B4 T- ?- U% n- p. g
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); . X, }1 C7 I\\" K, s; \
   
8. t=clock(),
   
9. / W. Q, O+ ^4 F5 K  Q4 g\\" g+ \$ p
10. oo{6 J$ O: s( w- d6 G
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. , ^6 @! y9 F) A8 r
14.     arrayfun[HFor("f"),x,y].Sum[]( `, v; H1 d* K# V7 [' |* y- |2 ~
    
16. };
    
17. # r  |3 g; A, g, D9 U6 b( {
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">4 J! s: Q1 ^& M6 h</font><br /> mvar:<font class=\\"jammer\\">% [4 B4 T- ?- U% n- p. g</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">. X, }1 C7 I\\" K, s; \</font><br /> t=clock(),<br /> <span style=\\"display:none\\">/ W. Q, O+ ^4 F5 K  Q4 g\\" g+ \$ p</span> oo{<font class=\\"jammer\\">6 J$ O: s( w- d6 G</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">, ^6 @! y9 F) A8 r</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">( `, v; H1 d* K# V7 [' |* y- |2 ~</font><br /> };<br /> <span style=\\"display:none\\"># r  |3 g; A, g, D9 U6 b( {</span> [clock()-t]/1000;

结果：
/ Y+ \% n, z, K* x) s5 M1008606.64947441
: n% P0 C% o  {& x4 q( G. X8 l! s) [0.735  秒
$ L+ u* G1 s6 M& p
& }8 d6 g, C" E7 q可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
& i6 s: u  z) }
－－－－－－－－

; G) }6 y, B1 M7 B1 q$ m& G从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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