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

forcal

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

1. //用C++代码描述为：; I5 g& e8 F  P5 ?5 {& `8 l( t4 \
   
2. s=0.0; . W+ \8 {& w. O3 _8 W+ t+ W
   
3. for(x=0.0;x<=1.0;x=x+0.0011) # _2 ?% Z6 r0 I2 e7 ~3 N
   
4. {( k+ m4 J% @$ b
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   * R# Q* ]6 a5 x1 a
6.    {
     P' O* G# n7 Y) o5 f% d- y# ?
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));/ S9 _+ y) Q* X* q. X
   
8.    }8 O! c! L: E5 x4 G8 S2 f9 Q
   
9. }

复制代码

Matlab代码：

1. tic
   : Q6 s4 q# M, {! s
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, @' P. y$ V& M9 e# v
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))/ A2 ]: x- O% s8 K* @
   
4. toc
   7 ?0 _( u2 [9 n3 r/ H8 R\" F' x
5. : i& s0 P9 X5 L! r4 c$ W! ]/ F
   
6. s =  M) `- I# U7 X\" l7 n9 m
   
7. 7 w  d) {- y6 g
   
8.   1.0086e+006
   5 K2 O; E; Z9 h# v, w
9. 
   . U5 |0 F9 `/ }3 I0 B; J
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 3 B\\" y) }/ V; p+ j; c& v# B
4. mvar:
   
5. % C5 |, \  K* u4 t' K
6. t=clock(),
   
7. 1 `( Y\\" K' A! u% T/ a4 A0 E
8. oo{; K& D. f: W0 {2 x5 J! h, v9 r
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],3 J# q8 w: o! v\\" n, }\\" \2 ]
    
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. ; r3 A: ]2 T- Q\\" ]5 n- d
14. };5 S& S; z1 g5 {; a! a
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">3 B\\" y) }/ V; p+ j; c& v# B</span> mvar:<br /> <span style=\\"display:none\\">% C5 |, \  K* u4 t' K</span> t=clock(),<br /> <span style=\\"display:none\\">1 `( Y\\" K' A! u% T/ a4 A0 E</span> oo{<font class=\\"jammer\\">; K& D. f: W0 {2 x5 J! h, v9 r</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">3 J# q8 w: o! v\\" n, }\\" \2 ]</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\\">; r3 A: ]2 T- Q\\" ]5 n- d</span> };<font class=\\"jammer\\">5 S& S; z1 g5 {; a! a</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
1 w& l9 q( r( [; z0.641

' }2 Y2 R3 r9 E+ |) Y( ?  ^Forcal比Matlab稍慢些。

+ n# R/ o9 x4 U; G* ]7 T% t－－－－－－－－－－
3 @& B1 u& Q" j* o% \6 U% G1 ?
3 G9 f/ ]" V- s9 h1 m再看循环效率。

" q2 v, P: ~4 bMatlab代码：

1. tic
   - n: C8 c0 a8 N) S! h
2. s=0;
   + }  X! _6 U; n2 v& r' W' {% C7 g4 Y
3. x = 0;
   * C  V\" y0 w( @
4. y = 1;- f+ ?  g* D; b4 f0 R5 N# p4 c
   
5. while x<1    9 C! u+ |: j\" w
   
6.     while y<2;        * r, p1 ~\" c' s9 K: a
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));\" i* f+ m/ |, \. w0 b# g/ }
   
8.         y = y+0.0009;        & e6 P& }9 [\" D: i
   
9.     end  L, A- j, v: v& z) x5 H- b
   
10.     x = x+0.0011;
    6 q: `: b* p0 m/ Z5 p3 f' Z: _
11.     y = 1;7 p3 S: _2 K( ~. y! U7 g( y% P
    
12. end, K$ L\" C1 {# J) u/ I, u
    
13. s
    # n1 a) P* z  \6 U- f. K
14. toc0 ]8 H+ t2 T% L+ d; m7 z/ @
    
15. # T  g4 v* h0 U* M  N( T
    
16. s =
    5 X. G# I' Q' b* `2 ^
17. : l% u! N- g! N; I. Q& u
    
18.   1.0086e+006
    8 ~( |0 S6 e1 p+ b6 e2 x5 o* e
19. / O8 H9 Y! |- L: M/ W4 w8 t6 l
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:\" s( Z* o* p8 D* b
   
2. t=sys::clock();( R! O7 ]+ D& j; @  {' v
   
3. s=0,x=0, 3 g+ ~4 Y  t* Q- |* u- K
   
4. while{x<=1, //while循环算法； + z6 m. S' n* P
   
5.     y=1, ( y9 g' Z. O1 J8 ]
   
6.     while{y<=2, 4 ~, d; @1 d/ ~! i
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 3 k2 x1 o3 |$ ~- p! L  I! R
   
8.         y=y+0.0009
   \" D1 d* R  o: L5 @: ^9 X
9.     },
   3 w. a' |2 `$ {& x( c7 A& [3 T
10.     x=x+0.0011
    9 }9 l5 d; f# M0 j
11. }, ! a7 a; M0 e0 c. D3 y; M9 E
    
12. s;
    1 q8 f$ B3 m; j  S
13. [sys::clock()-t]/1000;

复制代码

结果：
5 C  Y5 t% c# z- ^1008606.64947441
1 D7 O% R( @0 E) R8 Q4 W6 ^$ X0 K0.734   //时间，秒
  P, e: n; N; R
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

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

1. mvar:, o7 r6 a1 O6 g6 z$ p
   
2. t=sys::clock();  [* N% }, J# W! J/ a& u1 E
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   , [, E% K' Y# y. h3 m; _7 i; r/ k$ X
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   8 @# q. b2 p\" k
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0 A$ A- Q0 m3 _7 b! t3 e0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
5 ^* k4 Q- u+ I' {/ p. G
! X# M8 w  J( T% E  t2 {matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 E- V$ V( I$ E* I
2. tic, u4 Y& F) d% J# f5 [
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   4 x3 P& N, R, ]# K+ b
4. sum(sum(arrayfun(f,x,y)))
   8 W8 K- Y4 D1 B+ e  l\" B8 [( U\" I
5. toc; O) K& T# \' B( l8 `+ g8 `2 O4 \3 s
   
6. 4 ~  v. ~7 K  t8 M
   
7. ans =8 O\" i- L$ T3 E+ P* C
   
8. 
   4 M# g; g3 H, G: h+ a2 e
9.   1.0086e+0066 l+ P. K' c% L! H& a* w! T
   
10. 0 W/ y3 [( a. L' V) L0 R. z. O/ i
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3.   {3 }5 C\\" m+ U7 P& A
4. mvar:
   
5. , l\\" O* S6 Y( g1 }
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 4 _\\" m\\" I- N5 d4 |$ [# q
   
8. t=clock(),' n4 ^* W* Z# W\\" R8 q) O+ Y
   
10. oo{; X\\" z\\" b, }( U+ s/ b! o3 h. a
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( y5 Z4 S& s! n7 K6 f  l9 e
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 1 k/ ]5 N7 d0 f& ]
16. };
    
17. / Z0 k; I, p, y, M/ A' I\\" `' U) Z
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">  {3 }5 C\\" m+ U7 P& A</span> mvar:<br /> <span style=\\"display:none\\">, l\\" O* S6 Y( g1 }</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">4 _\\" m\\" I- N5 d4 |$ [# q</font><br /> t=clock(),<font class=\\"jammer\\">' n4 ^* W* Z# W\\" R8 q) O+ Y</font><br /> oo{<font class=\\"jammer\\">; X\\" z\\" b, }( U+ s/ b! o3 h. a</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( y5 Z4 S& s! n7 K6 f  l9 e</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">1 k/ ]5 N7 d0 f& ]</span> };<br /> <span style=\\"display:none\\">/ Z0 k; I, p, y, M/ A' I\\" `' U) Z</span> [clock()-t]/1000;

结果：
' N7 X' K  R, }, Z4 P5 u1008606.64947441
0.735  秒
; b/ j' ~- x& G& U) H
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
( H3 \, [, g2 s3 Q( w  F
9 y% K( G4 A% P, M+ Z! _9 ?－－－－－－－－
+ s/ U! `* Z1 V$ ]( S- b
6 x! w) u6 R$ [) M从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。