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

forcal

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

1. //用C++代码描述为：
   ; z7 K$ A\" i! h+ u* Q
2. s=0.0; ) g0 f# d0 U' E1 c+ _2 z) r- N
   
3. for(x=0.0;x<=1.0;x=x+0.0011) 3 E& U+ }. g6 _6 d+ k& a2 P2 O# }
   
4. {\" o; v( c/ E& o1 S3 \
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)/ h' D  B8 O6 j% I7 b, h6 H  X
   
6.    {! G: h- |3 N2 k* G' g  e
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   $ |* S, t8 K; h! U
8.    }
   2 q, `; t- v+ ?/ G% \8 M6 |/ _: s5 j
9. }

复制代码

Matlab代码：

1. tic\" W& K5 z' ]; q\" J, s\" V# {
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   & e7 Y/ ~( b, S  j
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   + a\" X4 C3 p  c# \/ d) i4 V5 {
4. toc. ~2 Z) r2 k. e% ^! Y$ X% s! ~1 c6 F
   
5. 
   + Y) C& B! q) h; e9 S& y. u! x- N
6. s =3 K/ z9 B1 x- A' a3 D5 M8 `7 i
   
7.   L% b8 P+ V& m2 o* X/ m8 S  y
   
8.   1.0086e+006. ?' m, d5 H. C$ Z/ n; E
   
9. 4 L4 U# \( O' L) G! a: ]0 Z
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];) a6 D& P3 k3 e, h
   
4. mvar:9 S+ X1 T  y' q\\" t2 {$ K# R  Q
   
6. t=clock(),
   
7. % m# {9 g6 V  N4 V: B* S
8. oo{
   
9. 2 U4 `; e6 ]/ v5 G7 x% H: i
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],! b& D) G) C6 A3 g; N% I
    
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. 6 x' J6 f\\" ]/ E- g! i- p7 U
14. };2 l: G% a+ O0 |0 v6 U
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">) a6 D& P3 k3 e, h</font><br /> mvar:<font class=\\"jammer\\">9 S+ X1 T  y' q\\" t2 {$ K# R  Q</font><br /> t=clock(),<br /> <span style=\\"display:none\\">% m# {9 g6 V  N4 V: B* S</span> oo{<br /> <span style=\\"display:none\\">2 U4 `; e6 ]/ v5 G7 x% H: i</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">! b& D) G) C6 A3 g; N% I</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\\">6 x' J6 f\\" ]/ E- g! i- p7 U</span> };<font class=\\"jammer\\">2 l: G% a+ O0 |0 v6 U</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
$ l. _! D9 O7 z2 z* ~3 h7 R& ?0.641

Forcal比Matlab稍慢些。
7 C% A+ \  B8 W+ b7 R7 @+ _
－－－－－－－－－－
7 k: ^7 Z* }; u4 f$ @
1 s. n- b. m$ [6 e( ~+ d/ F再看循环效率。
* d4 n! D# j1 k, u. \
Matlab代码：

1. tic
   7 {1 I# |3 \8 K- G* L
2. s=0;
   4 j. R6 a) K* J5 c/ `  c2 P8 v$ ^, y& k
3. x = 0;# D& U, y4 ^  M& L4 m8 Z
   
4. y = 1;
   \" z; D2 P7 R8 |% j9 Q$ D
5. while x<1    8 |; m- z' T5 w0 l  h3 m( h3 U. D
   
6.     while y<2;        
   9 |7 |6 r2 A& Z7 C
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 h* ?  J\" B3 F. X
8.         y = y+0.0009;        
   + h' r# v& V/ M7 e( l# m' f/ |
9.     end8 d, H* V$ m, m9 o; J) j7 Q
   
10.     x = x+0.0011;
    - a6 `- M1 Q. O# V
11.     y = 1;
    , P& N; s$ V+ V9 l1 _4 N8 i/ y, E) p3 \
12. end
    # J) E7 ^+ O. ^( ]. l. j
13. s
    6 k2 v+ i! b$ |- H- R\" R+ s
14. toc# p+ e$ S- u! X  o& {( L: q
    
15. 2 t5 L7 f6 q\" |) i! C
    
16. s =
      M+ M7 E' d) L2 d6 f& P! H
17. ; J$ `/ i  m8 u7 S5 ]
    
18.   1.0086e+006
    8 {! s\" Z* a+ h. T, ^
19. 
    9 q& }- F1 }( R
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   + n6 z% @7 g( L' ]; ~6 F4 H
2. t=sys::clock();
   % m# K+ U2 D. l- w) b, k
3. s=0,x=0,
   # s3 q7 T) R9 H+ n% @
4. while{x<=1, //while循环算法；
   . F! W6 s, i- d) Z3 W\" G
5.     y=1,
   ; N- r9 v7 A& p# b' X* R( S- i3 c
6.     while{y<=2,
   \" G5 W6 p7 t8 O6 r; D9 d
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   - M' Z2 I. j$ f
8.         y=y+0.0009
   1 z1 g2 z0 X2 d# v. c
9.     },
   - o1 D\" h: h+ A; u8 z7 W0 K; s; d
10.     x=x+0.0011
    1 i2 k4 n5 O7 U
11. },
    % H  u: L0 \7 ]' F0 a& K6 O- l
12. s;/ I. g* ?% P8 \; y$ W
    
13. [sys::clock()-t]/1000;

复制代码

结果：
8 e6 \* M7 |) X1008606.64947441
7 F  Z9 a5 i- K* |% _+ C( n; A) Z0.734   //时间，秒

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

6 ?1 T2 x( d' r7 H－－－－－－－

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

1. mvar:
   ! K3 z3 j8 R/ b1 Q
2. t=sys::clock();( |; w) N% ^' ]0 w  i3 b3 S
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % W9 q, C* E* R# M' `+ S) U
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ) b8 F) V2 G4 `9 t  A6 L; l
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
5 H9 }8 J- K. S
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - j, r4 u0 x/ a) O) u3 Z
2. tic
   8 i5 ^! c# ]) w+ Q\" x! m8 ^
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   - j- `# r/ c6 i0 `8 H8 e& u
4. sum(sum(arrayfun(f,x,y))); H5 e/ m5 n; J0 V' z, ]) ^
   
5. toc+ I2 }: g$ K7 g, Z
   
6. - f8 h  S% z; G( K' e
   
7. ans =
   4 A' W, J( D8 l& I. m7 o
8. 
   5 d\" e7 a2 D0 \' d
9.   1.0086e+006' r7 V1 o' U; l) U5 n
   
10. 
    ! g) z: c! j+ m7 X  H2 r4 z
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. . I! X* p% K1 _* ~3 l
4. mvar:2 g9 X& M1 P  ?# q) u! c
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 1 |& }# h+ v8 M+ l  b# X
   
8. t=clock(),* J, P0 U\\" B' O\\" S  Y  P+ X; J2 F
   
10. oo{) K* i- K- n3 a\\" x
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 5 H/ p, w1 h# Z3 e
14.     arrayfun[HFor("f"),x,y].Sum[]2 R, {$ w3 h) m! Q+ I, z2 F
    
16. };7 H9 s. Z/ r8 y% ^( }( ~3 L
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">. I! X* p% K1 _* ~3 l</span> mvar:<font class=\\"jammer\\">2 g9 X& M1 P  ?# q) u! 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\\">1 |& }# h+ v8 M+ l  b# X</font><br /> t=clock(),<font class=\\"jammer\\">* J, P0 U\\" B' O\\" S  Y  P+ X; J2 F</font><br /> oo{<font class=\\"jammer\\">) K* i- K- n3 a\\" x</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">5 H/ p, w1 h# Z3 e</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">2 R, {$ w3 h) m! Q+ I, z2 F</font><br /> };<font class=\\"jammer\\">7 H9 s. Z/ r8 y% ^( }( ~3 L</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒

3 a; F# m9 F1 X8 w6 ]可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

2 U9 W) H0 V7 Y, I. `1 _从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。