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

forcal

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

1. //用C++代码描述为：
   $ q& E. i9 y4 q3 I* ?( R
2. s=0.0;
   & @2 u3 p, C3 V' @5 d$ l; M8 C
3. for(x=0.0;x<=1.0;x=x+0.0011)
   1 J! ~6 i! S4 B0 o4 [, v, K
4. {
   + w/ z1 q* `\" L+ F% j& {: W: y/ U! o
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   9 a1 {2 ^/ @- \, W8 [
6.    {( D( u% `! w; B9 v9 c% N
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));. P9 g3 P. _! i5 c
   
8.    }
   2 f1 w( u/ p2 N5 J$ C0 H, o5 b; d
9. }

复制代码

Matlab代码：

1. tic
   ( H5 S\" ?1 r, o  x6 [9 q- l
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   . L4 R5 o+ F$ L* h9 X; b5 _
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))% G! n) h, v; S- g
   
4. toc
   5 X* I+ V2 S% Q+ Q. V. ^
5. 
   7 U- j+ @. F+ R) o3 ]6 L
6. s =
   7 T& j8 Y# @$ J0 N( D; A6 L
7. 
   6 w) _/ z1 A! p. p0 A( v4 I! {0 S% s
8.   1.0086e+0068 s! Y; F: o, b$ z) R4 q5 ~
   
9. 
   8 b, e/ U7 Q% g+ d
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];5 b\\" P# I$ d- d' t2 x! [8 N- Y# b
   
4. mvar:
   
5. 6 C8 k3 G9 p# Y4 D# ?
6. t=clock(),
   
7. . q5 B7 t2 W* O% d& i
8. oo{$ @1 \, A: x\\" y\\" A4 R8 g2 ?, S
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 5 z- v% |3 M% M- J$ o
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. 7 z# \\\" H7 g. d) l3 W( r+ a
14. };* q' g# r! A, f5 X
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">5 b\\" P# I$ d- d' t2 x! [8 N- Y# b</font><br /> mvar:<br /> <span style=\\"display:none\\">6 C8 k3 G9 p# Y4 D# ?</span> t=clock(),<br /> <span style=\\"display:none\\">. q5 B7 t2 W* O% d& i</span> oo{<font class=\\"jammer\\">$ @1 \, A: x\\" y\\" A4 R8 g2 ?, S</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">5 z- v% |3 M% M- J$ o</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\\">7 z# \\\" H7 g. d) l3 W( r+ a</span> };<font class=\\"jammer\\">* q' g# r! A, f5 X</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.641
) C: I2 M1 p8 [( X! W+ A4 ~5 _' E3 x
Forcal比Matlab稍慢些。
' o/ K! Q$ {3 O- G) ~
& ^4 D0 K) w; p+ {! |7 D－－－－－－－－－－

# S4 a- D9 V/ p1 h3 U0 w再看循环效率。
7 q$ [  |* v, ?
Matlab代码：

1. tic) T\" Z$ \+ }6 e3 E\" A
   
2. s=0;0 b8 w\" [$ ?  J, x' ?
   
3. x = 0;; X6 o: f$ W\" `. [
   
4. y = 1;
   \" r7 @& R9 U- b( ]  g0 s
5. while x<1    8 T) [5 J' {5 o6 l, j6 i
   
6.     while y<2;        * h+ T\" s+ t0 o3 I
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));- l- J5 z% \  i1 V
   
8.         y = y+0.0009;        
   $ t1 C/ h8 }( D3 V
9.     end( T7 I& K1 }, l\" W
   
10.     x = x+0.0011;, I\" P: k9 y6 S$ H4 _
    
11.     y = 1;  {- i* K0 _$ p+ d% H
    
12. end
    5 C2 H$ V' P\" ^; j6 T7 W7 A: j  O0 W
13. s
      b8 q! z, X% B
14. toc
    % a6 O6 B, |8 ?$ H8 z/ H% g
15. , y' h- X/ r7 B' [
    
16. s =% X8 p$ `6 a1 `% j3 Y: Y$ a
    
17. 6 A) \! j# y/ x# i; P
    
18.   1.0086e+006) Q\" J% M1 r* V2 f# f, R; ^\" x
    
19. & T0 N' O4 v8 Q0 Q) \7 Z
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   7 L% P7 v0 o! y+ R! B( n) A
2. t=sys::clock();8 w, D7 F8 i' N( c1 Z
   
3. s=0,x=0, 2 d7 a0 o1 D\" K
   
4. while{x<=1, //while循环算法；
   4 Z/ {\" g+ F8 {4 U! S
5.     y=1,
   1 Y9 r* M  y( z7 y2 r
6.     while{y<=2, . `! P) c+ L7 m5 L4 ~
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   , g$ F  x& |3 R; \5 W
8.         y=y+0.0009
   . r6 f0 U* R, {& T
9.     },
   2 D7 E) c7 \- H1 I
10.     x=x+0.0011
    8 ~7 y9 y, b5 G- H1 W4 Y& w
11. },
    3 q- a* N8 v$ B2 ^. X/ n3 I7 N; r) v) w
12. s;6 f; K% j! X2 @. s
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

5 F- m0 t' K# m* j2 }. |. p3 c$ I我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

8 ]0 J/ E; ^) k8 ]/ O" f2 o; M－－－－－－－

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

1. mvar:' O% ]; X; p$ A  f+ i; e
   
2. t=sys::clock();
   # N- s' o) ^5 m0 h4 T& ~- E5 |
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 h0 Q3 I; j! o3 {* y
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   3 g2 W\" @4 w4 P) a( n2 x! v7 W
5. [sys::clock()-t]/1000;

复制代码

结果：
" a" E4 |& M: G. p1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
1 ?) _8 C3 z, i( N% L4 y; Y
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));, z2 u\" w, t4 o* G, W
   
2. tic
   - E9 h1 h+ D! W/ c
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 n  s9 l( T* d- A+ u5 J) h7 Q6 }
   
4. sum(sum(arrayfun(f,x,y)))6 p* C6 s; ?4 b6 c$ K
   
5. toc
   & |2 B- a+ L- k3 b$ {: i! J
6. 
   1 x7 n* r6 Z$ ?% s( T9 i: `& e7 P5 F
7. ans =' x* J) g! _3 Q) B. a. R+ a
   
8. 
   0 B6 I2 E- J% U4 w7 y& e0 w/ N
9.   1.0086e+006& |% l9 g9 M* [9 t* O3 |
   
10. 
    4 p' q7 z8 c; r/ T$ j2 }9 w
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];) x1 ]) {. O6 f2 i4 S0 t8 Y; j3 z
   
4. mvar:\\" s8 k- |\\" A5 m; t# B. Q\\" z
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 4 H; z) F\\" o' [. Y. J
   
8. t=clock(),- V( W* U7 Y3 r4 l
   
10. oo{3 ?' \9 ^' T2 f% y
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],9 K# @8 r; ]7 ^1 B, K3 ~0 q8 g
    
14.     arrayfun[HFor("f"),x,y].Sum[]. N, d& j+ ^. Y) w\\" M  L
    
16. };6 [, ~% \0 e& C\\" C; D+ i6 W0 y6 R( ]
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">) x1 ]) {. O6 f2 i4 S0 t8 Y; j3 z</font><br /> mvar:<font class=\\"jammer\\">\\" s8 k- |\\" A5 m; t# B. Q\\" z</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">4 H; z) F\\" o' [. Y. J</font><br /> t=clock(),<font class=\\"jammer\\">- V( W* U7 Y3 r4 l</font><br /> oo{<font class=\\"jammer\\">3 ?' \9 ^' T2 f% y</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">9 K# @8 r; ]7 ^1 B, K3 ~0 q8 g</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">. N, d& j+ ^. Y) w\\" M  L</font><br /> };<font class=\\"jammer\\">6 [, ~% \0 e& C\\" C; D+ i6 W0 y6 R( ]</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
4 A( P1 F; T0 t" {& O* u& w0.735  秒
. S( v% p; y5 P, y: F0 l) R
9 C  u2 y! |' l' e/ y* q可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

2 K) G: Y+ s' w. D& R7 @& y5 w) T－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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