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

forcal

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

1. //用C++代码描述为：
   + b* O1 a3 C, k' m& [6 H$ h- l
2. s=0.0; ( X! G4 a! O& {# Q, P( R6 [9 B, \
   
3. for(x=0.0;x<=1.0;x=x+0.0011) ! k' }0 o3 d. K! y
   
4. {# k. E\" |( f; v) }, I
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)3 w& @: o: k# P: [, e
   
6.    {
   3 A! L& b9 I+ K0 B6 _5 |1 t4 H5 s
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));, r2 a; ^  H\" @! j4 p4 n8 W* _2 X
   
8.    }! g5 N% |/ Z& S9 s
   
9. }

复制代码

Matlab代码：

1. tic* O. N) ?: E. J
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ' I; k+ h- b1 H1 b* I\" q# Z
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   - U3 {4 a. B; q4 H\" C' L
4. toc
   ( A\" i0 A* T0 H; [9 ^! E1 F' Y
5. ' f7 u: w4 F$ f  s: a
   
6. s =
   : s& }8 f& v. i' S9 B# M
7. 
   : i3 I. w8 z1 _\" r6 Z
8.   1.0086e+0060 y  J7 Q3 ?$ p
   
9. ' U7 {8 G- A' }1 u
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];- X3 }6 p  d- s# f3 J- X0 l% J( R
   
4. mvar:
   
5. 5 }8 d# W& w+ e' ~
6. t=clock(),
   
7. 7 ]; `2 Y6 P! x% y
8. oo{
   
9. ! D9 e, y6 y: Z5 _3 X
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],; ]: M' L5 ^* `. I7 c* F* |- Y# U
    
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]( Q  ~: D% |& X' r, k7 O6 }
    
14. };( }0 R  m' |/ q) v. m9 y' A% f
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">- X3 }6 p  d- s# f3 J- X0 l% J( R</font><br /> mvar:<br /> <span style=\\"display:none\\">5 }8 d# W& w+ e' ~</span> t=clock(),<br /> <span style=\\"display:none\\">7 ]; `2 Y6 P! x% y</span> oo{<br /> <span style=\\"display:none\\">! D9 e, y6 y: Z5 _3 X</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">; ]: M' L5 ^* `. I7 c* F* |- Y# U</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]<font class=\\"jammer\\">( Q  ~: D% |& X' r, k7 O6 }</font><br /> };<font class=\\"jammer\\">( }0 R  m' |/ q) v. m9 y' A% f</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
6 \3 x8 ]& [! s0.641
* g' l# x- y8 |6 L* S0 I" T. X
Forcal比Matlab稍慢些。

$ z3 w4 D2 t* `9 D+ `( ?－－－－－－－－－－

再看循环效率。
. f) k# l. k$ _
Matlab代码：

1. tic
   9 ~# k3 A- S7 \7 O) v9 e
2. s=0;7 q0 `/ l& Y: ]& s% d! ]
   
3. x = 0;4 b% S, B( W7 i' ~
   
4. y = 1;
   2 R  ?5 b2 `7 ], J2 H2 j3 S
5. while x<1    ! |4 j+ @) ~  J
   
6.     while y<2;        8 q; Z- T/ E' O\" E/ E  e
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 ^\" Y3 l3 {, b+ v' @: _
8.         y = y+0.0009;        
   + I% i+ T' e$ d! q4 l' V9 B: e
9.     end
   ; V! d9 ]2 P; ^% V/ Q# M. G& _\" l
10.     x = x+0.0011;
    ( f1 v) K\" d4 m. |: V
11.     y = 1;5 Z: \$ I, g+ b/ t6 t, \5 b
    
12. end$ ^- W% I/ f( \* P1 r
    
13. s3 n8 U0 t! y\" S2 `
    
14. toc. U; W$ I6 [6 {\" H
    
15. & I6 A& ]! R: W4 O  t! ?
    
16. s =( \& [1 R7 P4 C+ U' I8 Y1 C1 F# M
    
17. 1 q1 R5 i. e\" J& H0 K
    
18.   1.0086e+006$ B2 G: Q2 L: o4 L0 Z% {\" s6 T+ S
    
19. 
    1 }& }8 K/ I8 b& J  C
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   ( I- y, X7 p& T& b' Q& {
2. t=sys::clock();3 M! `8 }+ ]7 r6 E- n
   
3. s=0,x=0,
     M( O9 w# S. [: E
4. while{x<=1, //while循环算法；
   : i$ D0 J  d, V- `1 c0 l& o8 c( Q
5.     y=1, ' x2 N$ {6 O' ^' r: r
   
6.     while{y<=2, : x4 N6 q6 U: j' G+ F3 x
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   , y+ }\" g1 P8 Q) {2 R
8.         y=y+0.0009
   \" Y) l- W% S+ G; a
9.     },
   4 U% {; V' I* u# u8 {\" c\" T
10.     x=x+0.0011 + ~8 l5 G\" I* ?! B( ^# K
    
11. },
    \" {2 m8 x) v% L\" `! t
12. s;& n6 y; f& K2 Y\" G4 x- f- _* H8 Z' F* L
    
13. [sys::clock()-t]/1000;

复制代码

结果：
% J5 p, O, p6 e! F- r% N9 y1008606.64947441
0.734   //时间，秒
+ D2 f' t% U( N! |) B
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
, n8 d( ?' m/ j
－－－－－－－
5 I" K# }8 o5 a2 I: Q# w* b" N
0 j- y! E; R# R8 X* VForcal中还有一个函数sum专门进行这种计算：

1. mvar:& n1 L/ |, j$ j5 g' T' p4 x
   
2. t=sys::clock();
   ' [5 s( B6 @& k5 _9 N6 S
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 5 X2 {+ U3 b9 ]) |% o# L- ]( t' ?
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   & u  g2 `9 g$ B. v6 Y2 [
5. [sys::clock()-t]/1000;

复制代码

结果：
4 w% i8 c  i: t* S( e1008606.64947441
0 I0 L' B% X: M; ^- z0 Z9 m6 s0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
% y+ k' j: v+ |& n2 ]: E/ ~2 E, w
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - g' N3 V  c- i2 g0 M/ N# l: x
2. tic! w  n, n9 {$ B5 ^\" p
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);/ H7 e8 |: f+ y: U, \
   
4. sum(sum(arrayfun(f,x,y)))
   : a  E4 D! C5 J; v, A) v
5. toc
   , q1 [# l0 ]5 z; m4 Y
6. & C\" V: i2 g- z1 K2 |1 x' U
   
7. ans =  [. @* X& C3 k- A: H
   
8. ! d3 n: E2 a2 F7 y/ m7 k& y
   
9.   1.0086e+006; \& v* e+ x5 t9 ?  g4 C; S
   
10. 3 F9 B& W5 }6 d
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];. W: A( p# C* s1 \# }8 q. Z& I
   
4. mvar:
   
5. ; \\\" B\\" m2 j, q% T2 I- }
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ! v6 ?5 v1 Z& O0 ?
   
8. t=clock(),% ]% d8 S: O. N% Y  Z
   
10. oo{* t2 W: P) P/ R) L3 Z
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. ! m, {: Z4 ?& l* P( N% `/ v+ n+ k' v
14.     arrayfun[HFor("f"),x,y].Sum[]) n/ Q1 d( N  E9 d6 P3 E
    
16. };
    
17. ( E4 z9 O5 m. n7 I
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">. W: A( p# C* s1 \# }8 q. Z& I</font><br /> mvar:<br /> <span style=\\"display:none\\">; \\\" B\\" m2 j, q% T2 I- }</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">! v6 ?5 v1 Z& O0 ?</font><br /> t=clock(),<font class=\\"jammer\\">% ]% d8 S: O. N% Y  Z</font><br /> oo{<font class=\\"jammer\\">* t2 W: P) P/ R) L3 Z</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">! m, {: Z4 ?& l* P( N% `/ v+ n+ k' v</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">) n/ Q1 d( N  E9 d6 P3 E</font><br /> };<br /> <span style=\\"display:none\\">( E4 z9 O5 m. n7 I</span> [clock()-t]/1000;

结果：
' p  u% n6 Q& T/ |" r8 z1008606.64947441
7 f! _' U. a7 N0 Z) \  R# R) F; X0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

$ M. A2 L4 t( w. R－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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