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

forcal

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

1. //用C++代码描述为：
   - Y8 s, s9 C\" N; ?
2. s=0.0;
   0 R1 W0 p1 @2 N; |3 Q6 ?1 n
3. for(x=0.0;x<=1.0;x=x+0.0011) 0 o) ?+ C8 R, U( A% }
   
4. {* j- l; j6 c% m\" |: K) K
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)3 ]0 l\" k% K7 \4 d  [# w; b) ]
   
6.    {1 c: h  _9 Q* u\" w% p\" [
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));3 ^) f% k2 ]: |* c$ P- P4 q
   
8.    }
   , `$ ?3 a; h) C, |6 D
9. }

复制代码

Matlab代码：

1. tic
   . n8 K+ \4 U! F+ A6 F  i; ]
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   1 s9 b- a' y9 z) R9 Z& Q
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   4 O6 S7 z# ~8 l0 M
4. toc7 i4 O9 [2 a; N+ O/ M8 b% G- T
   
5. 
   , X$ v- n) u- B! ?' J+ ^
6. s =
   ; {$ {1 g) B' ?+ |$ R7 p' o
7. 
   ) ?# L) u! c, }5 `# C! R
8.   1.0086e+006
   % c0 Z5 F8 C, u1 Y$ ~
9. 0 W9 p# o# m; m8 P
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 0 S0 a) W, q/ c7 ~3 k9 b8 A
4. mvar:
   
5. , e6 m) a! F: A9 W6 ?2 Z
6. t=clock(),( ?2 G3 u5 w1 v; ]9 k1 O
   
8. oo{\\" x  r9 Q, }. l; O
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],/ n2 ?  u( ?! [9 z6 N
    
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. 1 u2 f. T+ k$ x( G. j
14. };
    
15. 4 F\\" U0 ?$ j: k3 K8 O* Z* k+ S
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">0 S0 a) W, q/ c7 ~3 k9 b8 A</span> mvar:<br /> <span style=\\"display:none\\">, e6 m) a! F: A9 W6 ?2 Z</span> t=clock(),<font class=\\"jammer\\">( ?2 G3 u5 w1 v; ]9 k1 O</font><br /> oo{<font class=\\"jammer\\">\\" x  r9 Q, }. l; O</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">/ n2 ?  u( ?! [9 z6 N</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\\">1 u2 f. T+ k$ x( G. j</span> };<br /> <span style=\\"display:none\\">4 F\\" U0 ?$ j: k3 K8 O* Z* k+ S</span> [clock()-t]/1000;

结果：
9 `; {) P, L8 t8 a, `& N1008606.64947441
, L$ s7 N  x# X. R0.641

( A1 d7 H: k5 I8 s! K0 xForcal比Matlab稍慢些。
9 \( x. x9 H& G! s
/ ?- ], G; V/ G& N0 p0 G－－－－－－－－－－
7 }0 m1 ]% m# R7 B
再看循环效率。

Matlab代码：

1. tic2 y0 F* V5 N; H: H) `
   
2. s=0;
   % R6 _$ p+ M, g6 H
3. x = 0;
   : p& z. ^5 Z- g; R4 Q1 E
4. y = 1;\" Y4 Z6 c1 Q0 N1 `( `7 V
   
5. while x<1    9 f5 `% s% R5 H: Z2 L; d
   
6.     while y<2;        
   / P1 R2 K7 i& [3 b2 N' K7 k6 d
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( _# R* H: g. W* O
8.         y = y+0.0009;        
   9 f$ f5 U5 ]  y! O- F$ |
9.     end
   ( c+ L\" L; o7 C! B$ }+ {% e
10.     x = x+0.0011;1 Y1 Q7 I8 B. _' O' a
    
11.     y = 1;
    * K( S+ c, O1 g$ i\" A9 [; a
12. end9 M. `: z/ X9 a  n' n
    
13. s
    ( j/ `$ s$ d! n, @  I7 V; M' z' b
14. toc6 Q9 @( c* v! I' I
    
15. + M5 v. L; `' Y8 l% L
    
16. s =
    ! L! q( P5 n* i\" E4 L
17. 
    ; F* }0 c$ v1 b: t7 ]  e
18.   1.0086e+006
    & R\" F& I) g( ~1 L# {1 o+ w
19. 
    4 }* \\" s\" S6 J\" U7 g3 y. h
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:9 x8 {6 s/ j* L* K  y: [
   
2. t=sys::clock();
   ' ?* a- K3 [+ V$ L; L7 Z8 ~* S9 ~
3. s=0,x=0,
   + b9 m2 R% I9 T; A) Z5 F& T
4. while{x<=1, //while循环算法；
   9 Y2 w& f% i5 A4 g7 n2 Z. W3 R
5.     y=1,
   & N\" [8 e$ P  J1 i/ o# H
6.     while{y<=2, ) D9 v$ S  W- C- r( b& N
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   / w& y/ y) |, c( u
8.         y=y+0.0009 & a7 f- Q7 a' y2 }* ?2 d
   
9.     }, 0 i! N5 |% R7 f0 C
   
10.     x=x+0.0011 & \* @) p- r( M. W( e, c% w, O
    
11. },   N% P& U$ g8 h9 D; ~) x; k\" i
    
12. s;
    2 s7 M0 U: o  G/ U2 x
13. [sys::clock()-t]/1000;

复制代码

结果：
2 d  R4 R  l4 k% U/ p1008606.64947441
0.734   //时间，秒

, T# |7 X1 n* @( r( q% E& q$ W+ f我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
$ |/ E3 C: O. p! x6 }
$ c+ p1 H' I6 S# c' f- p1 c9 o9 q－－－－－－－
2 v* Y9 c5 N$ @5 s$ Z/ ]8 c: u  Z
  W: P# ~& }5 e. G- }0 V, rForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   6 H( N+ E3 P' F- k+ D
2. t=sys::clock();% I. y\" W6 Q9 F8 H
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 s0 ~0 s\" @) Y! Z\" J
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   0 m3 Q2 j* _1 i) z2 s$ J8 C' V+ @
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

6 M$ ^8 `/ x; t& f9 ?; k" M6 `7 kmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));4 m1 ]/ E, `( s. F: K0 m$ t/ X
   
2. tic
   $ `0 Y/ l/ B- D+ a
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);9 J& c! u, j8 v& O8 q- X\" I
   
4. sum(sum(arrayfun(f,x,y)))* L4 i8 C% @* r% @0 V' t
   
5. toc) M: Z2 W' F0 J1 \4 v6 i
   
6. : R  |, [5 v\" X\" X
   
7. ans =2 j* V8 \8 [# W, f  G7 K; k# Q
   
8. 
   4 W3 A: ~/ h) S0 j: G
9.   1.0086e+0067 E7 G8 q: C  h3 v# |. ~3 z
   
10. 
    - w: L/ _' \$ B$ o3 B
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];1 Y% p3 d; y' i! ^0 M, m3 {) \
   
4. mvar:8 N8 |* U% F6 |3 H* h
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); \\" R2 ^  w1 u0 T9 {  T
   
8. t=clock(),
   
9. , b6 L( m$ ~* N* x
10. oo{
    
11. \\" m1 ?1 Y, {3 G- Z- p
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. % g6 L7 q2 {$ H& M* r
14.     arrayfun[HFor("f"),x,y].Sum[]% }& ^% h& K0 P+ U9 C
    
16. };2 Q\\" p' Y- }) P; W7 n; k: t3 T
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">1 Y% p3 d; y' i! ^0 M, m3 {) \</font><br /> mvar:<font class=\\"jammer\\">8 N8 |* U% F6 |3 H* h</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">\\" R2 ^  w1 u0 T9 {  T</font><br /> t=clock(),<br /> <span style=\\"display:none\\">, b6 L( m$ ~* N* x</span> oo{<br /> <span style=\\"display:none\\">\\" m1 ?1 Y, {3 G- Z- p</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">% g6 L7 q2 {$ H& M* r</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">% }& ^% h& K0 P+ U9 C</font><br /> };<font class=\\"jammer\\">2 Q\\" p' Y- }) P; W7 n; k: t3 T</font><br /> [clock()-t]/1000;

结果：
6 e& T/ }. `8 }3 M1 B9 Y1008606.64947441
0.735  秒
5 D) @1 I- |+ i% C7 Q% w
% E6 H% ]2 d3 j* ]- n6 n* l可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
. }7 k3 I$ S9 H+ b* m3 c
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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