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

forcal

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

1. //用C++代码描述为：
   7 L8 g6 w, L/ S
2. s=0.0;
   + W& g/ _; B: ~4 z6 S6 O
3. for(x=0.0;x<=1.0;x=x+0.0011) ( _4 \& i* u9 ~, w* J  M1 ?
   
4. {! d% O\" y' U\" t& v$ N! M
   
5.    for(y=1.0;y<=2.0;y=y+0.0009), {& N3 O( g3 ?, D
   
6.    {
   8 O. r1 y( R& S  C% }. ~
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));1 P0 t; C\" l9 Z1 p- o' r
   
8.    }9 Q3 y$ B3 B. q$ {\" B. a; U; q
   
9. }

复制代码

Matlab代码：

1. tic
   8 d3 s: o+ S  T( N. a8 A
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
     |\" l8 H  Q& m/ ^  Q( K. N( n
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ' y8 |! x7 U2 \! S( e5 H
4. toc
   * R) n1 l  x* C0 v8 n, Z0 v\" C
5. 
     t; C3 S6 o9 X( Z# A+ J0 z* s3 R
6. s =
   + ]9 q0 j' }7 l7 }' O8 V
7. 
   0 |1 J* k) C+ T  {( ]  M
8.   1.0086e+006# t6 e; u# c4 K
   
9. 
   2 W1 m+ v1 o* m% u% u& D5 H& o$ l
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];1 h+ X% t\\" |( l3 x$ U, [; v
   
4. mvar:
   
5. $ V7 g  A2 y8 c2 F7 m/ b! i
6. t=clock(),+ L. l, v5 ^) `, r
   
8. oo{
   
9. . `( \9 A7 V. g# c
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 6 P( n4 `( H; R
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. , W& ~9 i9 ~, h- b/ {
14. };
    
15. : ~: c$ n( \& s8 X8 y- I9 R  t3 ?( Y
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">1 h+ X% t\\" |( l3 x$ U, [; v</font><br /> mvar:<br /> <span style=\\"display:none\\">$ V7 g  A2 y8 c2 F7 m/ b! i</span> t=clock(),<font class=\\"jammer\\">+ L. l, v5 ^) `, r</font><br /> oo{<br /> <span style=\\"display:none\\">. `( \9 A7 V. g# c</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">6 P( n4 `( H; R</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\\">, W& ~9 i9 ~, h- b/ {</span> };<br /> <span style=\\"display:none\\">: ~: c$ n( \& s8 X8 y- I9 R  t3 ?( Y</span> [clock()-t]/1000;

结果：
3 r# O2 l, w2 T; B; s1008606.64947441
  _3 b. Q6 w/ ]0.641
& U* o5 Q5 j; o- y
Forcal比Matlab稍慢些。

8 Y( u1 E3 B2 H/ F－－－－－－－－－－

9 q/ B4 E% F! k* ^2 u& f再看循环效率。

Matlab代码：

1. tic
   ) w  a  N9 I/ l( H6 O' X
2. s=0;
   7 [3 @! H4 g7 R6 S0 e7 N: h+ N
3. x = 0;# `$ I: E6 @  y! E$ ^: y
   
4. y = 1;
   & V7 \$ g1 j4 Q5 s& i
5. while x<1   
   . V) V* N1 c4 a0 f; t2 y$ {( y+ d
6.     while y<2;        # q: ?6 V. ]\" Y
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / a, `! i\" i! _( J
8.         y = y+0.0009;        
   7 A  V+ V! S7 s7 _
9.     end
   & T9 y# b' Q, W4 L
10.     x = x+0.0011;2 m\" L9 Y& F: }, j; F
    
11.     y = 1;, a) g( k( m, a5 @# p& E& L
    
12. end( e$ ~+ o/ C; \: N9 A2 `) T
    
13. s\" L& K  U7 S3 n' C7 D8 J
    
14. toc. K( R; }4 D4 o3 |% n3 z/ D) q
    
15. 
    4 I* v3 t9 Z% V$ `$ j+ V
16. s =) M\" P- W) {) K1 F\" ~\" j
    
17. 
    8 _( I  b! z! d6 q3 [2 {
18.   1.0086e+006/ {9 u. ~# g) [- c: Y
    
19. ; }) B7 c0 f  `; w  |
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:  A& }. t) H, a  G
   
2. t=sys::clock();
   ! y* m8 u9 q7 V
3. s=0,x=0,
   . E3 I! J; i. L! P% j7 Y
4. while{x<=1, //while循环算法；
   , y- P0 w; p5 g
5.     y=1,
   \" G6 {& a) q% S
6.     while{y<=2, 3 W& S! v$ C. d; H6 n+ \
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   . H' y+ I; ?2 {  J% U
8.         y=y+0.0009
   , W6 ~- H8 R2 _
9.     },
   - |# A/ B! v- u( q: j- }9 k0 r3 \
10.     x=x+0.0011 9 b+ U  U+ ?3 z
    
11. },
    & s8 n& Y+ m, Z. O; q* B1 H
12. s;
      F5 q4 b& E' i1 w( Y) ^5 E  z
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
& r# N; w. F. d0.734   //时间，秒
' ?5 _# p/ y4 l9 i; }, w
/ {9 L. M( j( Q0 u我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
3 O3 P5 q* O6 r$ H* A9 G
－－－－－－－

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

1. mvar:$ q6 _3 F\" q2 a0 I4 B. x
   
2. t=sys::clock();. }) M% ^) v% Z1 p/ l7 x
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( P\" \8 Z# a: p8 Y, B% I5 X
4. sum["f",0,1,0.0011 : 1,2,0.0009];' j6 t: e4 W$ {' F( U5 t
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

9 D' m8 h- g& R7 f2 S) M# cmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) I  _  z) u\" b' R\" k8 U
   
2. tic9 g\" @) _2 z8 `# C4 |/ S7 j$ o
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   % q5 R\" O# {: i9 H. o% F\" i+ g
4. sum(sum(arrayfun(f,x,y)))
   + i. r: `0 E+ g, l- ~
5. toc
   - h7 k2 J* d. T5 B' x
6. 1 L  P  {+ r! v7 ?
   
7. ans =3 V% ~$ I. n) R+ q4 D
   
8. ( ~1 p6 K% c6 B
   
9.   1.0086e+006) e& }3 V$ D9 n
   
10. 
    1 }7 Z, n% @* Y$ m9 c
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];7 N* L6 _6 a) \. }. Q
   
4. mvar:9 r/ S5 s# y& B, U% d3 K
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); $ L: t# F- i* h$ _+ s+ W# D9 j
   
8. t=clock(),6 i& ^4 {3 U\\" t+ F
   
10. oo{
    
11. 6 K1 b+ z0 f# F$ c
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. / P' z7 o4 Q# O
14.     arrayfun[HFor("f"),x,y].Sum[]' o3 ^( d# \: c3 M( S
    
16. };2 a7 X% `3 X( K  {
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">7 N* L6 _6 a) \. }. Q</font><br /> mvar:<font class=\\"jammer\\">9 r/ S5 s# y& B, U% d3 K</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">$ L: t# F- i* h$ _+ s+ W# D9 j</font><br /> t=clock(),<font class=\\"jammer\\">6 i& ^4 {3 U\\" t+ F</font><br /> oo{<br /> <span style=\\"display:none\\">6 K1 b+ z0 f# F$ c</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">/ P' z7 o4 Q# O</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">' o3 ^( d# \: c3 M( S</font><br /> };<font class=\\"jammer\\">2 a7 X% `3 X( K  {</font><br /> [clock()-t]/1000;

结果：
0 O$ ^1 {- j% y$ [+ d( h4 R1008606.64947441
( j, `% X8 r3 J7 {+ w0.735  秒
+ l+ S8 `% A' ?; ]9 V' u
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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