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

forcal

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

1. //用C++代码描述为：
   1 x4 m% l$ {1 s: h: e) m/ y
2. s=0.0;
   ; D9 l8 D& c1 d/ }* q
3. for(x=0.0;x<=1.0;x=x+0.0011) 1 P5 m/ k5 x9 |& i% }0 N3 H6 f& z1 S
   
4. {1 ~1 J1 {) \  J5 s5 e8 v7 d# b
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)2 P. m1 o* z* J
   
6.    {. n; K! A0 f+ I* x- y
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));; k* M; h; O) D) n  T5 [
   
8.    }
   # x9 O7 B9 M3 q
9. }

复制代码

Matlab代码：

1. tic+ E/ {8 T8 v7 |
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ) ]\" @% H- ]+ l\" V6 d
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   % L& t  X# H, N
4. toc& f5 v7 v7 Z4 e3 R% M
   
5. 
   + Z8 [$ \4 H6 W  g+ }6 l
6. s =& w8 Q: f\" g. m+ E: l* g# x
   
7. 
   ! o* i3 w# n: e. w9 o
8.   1.0086e+0069 g' \* A+ S; t  N$ k+ Z3 P
   
9. 
   : I5 w4 S, u  @
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; h# `9 E3 i; m1 R0 V
4. mvar:0 I: {0 d7 j3 z' G5 a4 D$ h. U# `9 S- e
   
6. t=clock(),; @7 f$ f8 W7 L) A+ F( F4 V  V. |& @
   
8. oo{1 m* o0 u% O\\" e( B# g
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],2 g8 b. V; k! R: 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], B) A3 P0 L* |9 `( j
    
14. };: h7 e; q: K% O* j$ ~
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; h# `9 E3 i; m1 R0 V</span> mvar:<font class=\\"jammer\\">0 I: {0 d7 j3 z' G5 a4 D$ h. U# `9 S- e</font><br /> t=clock(),<font class=\\"jammer\\">; @7 f$ f8 W7 L) A+ F( F4 V  V. |& @</font><br /> oo{<font class=\\"jammer\\">1 m* o0 u% O\\" e( B# g</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">2 g8 b. V; k! R: h+ R</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\\">, B) A3 P0 L* |9 `( j</font><br /> };<font class=\\"jammer\\">: h7 e; q: K% O* j$ ~</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
% }2 {6 J5 R% X+ G0.641
; d% p5 |* Q" B7 Q4 f
Forcal比Matlab稍慢些。

－－－－－－－－－－
$ O5 ?5 n3 F+ `$ |5 e! k- U6 w8 V
5 Q  R. }2 U$ E( o5 P; N再看循环效率。

Matlab代码：

1. tic
   5 ]& w0 j1 d' c: D
2. s=0;5 a2 o! V: [4 a2 M
   
3. x = 0;6 {. g9 w7 ^\" [' }0 _: r\" ~
   
4. y = 1;0 {- L5 K1 o: X* m
   
5. while x<1    6 w' o! q* v) g: I; d\" t
   
6.     while y<2;        
   8 C$ B% ~5 O* X, \& E
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   : i8 j% P. {! [  d2 \\" r0 V9 P$ R
8.         y = y+0.0009;        
   ' N; `3 I! T9 R8 h2 N% p4 F- z! `
9.     end9 a  b2 `1 Y7 n\" M
   
10.     x = x+0.0011;1 i( L& J' d: G1 C
    
11.     y = 1;* D) v9 S& k! o; d4 R6 ]5 Y- w! A
    
12. end1 g% d% X1 ]' ]2 U4 c6 m
    
13. s' B. X) r) L  _  K6 M, `
    
14. toc# G9 q( d1 m: Z: l: r. `8 ~, N
    
15. 
    * `/ d, E/ x/ F1 U- ?8 `3 P, N( l
16. s =
    ; Y; l! B5 j5 x( o
17. 
    1 e/ o6 w9 O, I& G
18.   1.0086e+0062 Z( J, W4 B5 |\" A6 M; K\" J
    
19. , b7 Q# P, n: a. |+ N% K\" P' R3 J
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   & O3 M2 g$ V5 k1 i) _. t# |
2. t=sys::clock();
     [9 `' }6 _' J9 h& d; z
3. s=0,x=0,
   ! `) {/ d2 c. M
4. while{x<=1, //while循环算法； ) r1 T# H- ?$ V0 a0 t
   
5.     y=1, ; G* c! ?, {4 S: G/ D
   
6.     while{y<=2,
   5 Y0 i  \3 i# q& N& B/ w6 V' ?
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   3 i; v2 `2 K/ W# q7 q7 p( B
8.         y=y+0.0009 ' F% O9 A: l+ O0 L/ o
   
9.     }, - R) d+ r5 \2 |* m9 X9 F/ U; t& T* n  D
   
10.     x=x+0.0011 - W' A6 z5 S  n0 E0 y+ P1 J& Z
    
11. },
    7 Z$ F- q, `: L. b1 {) k
12. s;; @' s: \6 Q' n4 B
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
5 M0 _: E! |+ c0 x( p9 N* R0.734   //时间，秒
) L, c' N. ~* R( e
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
5 F4 e8 l$ F2 N' C/ k4 Y8 r
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   8 u6 _: \) Z  ^- k6 x\" V( y, j6 F
2. t=sys::clock();  i5 s8 M# ~- i
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 9 D\" K, E' ?5 M
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ; m2 H' J6 o$ t! Z
5. [sys::clock()-t]/1000;

复制代码

结果：
7 Y# i! `& ~7 M/ M- v$ E1008606.64947441
0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));' [6 S8 k; }6 Y/ F
   
2. tic
   . y0 a9 F& j- g# g! X9 d2 Q9 _. L' w
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   $ Y- H& P2 B% _
4. sum(sum(arrayfun(f,x,y)))/ j& g2 ~+ |( p* N
   
5. toc
   ) x9 P7 I+ d3 j5 U) C
6. 
   , `! d9 w! I( j7 Z( Z
7. ans =
   & i\" c. c/ t; Q1 r8 H. e$ m
8. 2 ]6 w3 w$ q% Z4 C5 d0 {
   
9.   1.0086e+006- y3 C, E+ h9 K  c4 `5 ?- y
   
10. 
    : A8 i9 R' [) k4 Y% u9 J2 Z/ f
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ( c# n6 y! g\\" n9 G7 U
4. mvar:9 b; Q9 w* Q2 I6 s6 t/ t
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 8 P7 l& }6 d' ~
   
8. t=clock(),
   
9. ' ^/ m9 q8 D8 i& [, X
10. oo{
    
11. ! G+ L1 }8 N0 r2 e6 ^: K
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. $ S3 b5 A) g5 A7 Y- b
14.     arrayfun[HFor("f"),x,y].Sum[]: f8 {' l+ A$ M/ B0 ~( Y
    
16. };) Y5 j. ], q$ L; A; r0 w5 p
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">( c# n6 y! g\\" n9 G7 U</span> mvar:<font class=\\"jammer\\">9 b; Q9 w* Q2 I6 s6 t/ t</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">8 P7 l& }6 d' ~</font><br /> t=clock(),<br /> <span style=\\"display:none\\">' ^/ m9 q8 D8 i& [, X</span> oo{<br /> <span style=\\"display:none\\">! G+ L1 }8 N0 r2 e6 ^: K</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">$ S3 b5 A) g5 A7 Y- b</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">: f8 {' l+ A$ M/ B0 ~( Y</font><br /> };<font class=\\"jammer\\">) Y5 j. ], q$ L; A; r0 w5 p</font><br /> [clock()-t]/1000;

结果：
: v% S4 @1 c; M  {5 Q, L1 E1008606.64947441
( {% \+ H% f7 o, y3 t0.735  秒
6 e1 m+ X( V- A  f7 I
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
( q9 `9 L8 P' i9 q; s
- v: Q! b* B# T: D－－－－－－－－

# o$ P  X" u/ `9 p; v# \2 y从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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