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

forcal

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

1. //用C++代码描述为：+ Z% g( u2 L1 P( I& N
   
2. s=0.0; , u' \\" ~2 r& J; E9 l7 v
   
3. for(x=0.0;x<=1.0;x=x+0.0011) \" H* M, B. W7 v
   
4. {
   1 ]9 c1 ]  m6 Y* a
5.    for(y=1.0;y<=2.0;y=y+0.0009)* p5 X- ?3 |: O% U
   
6.    {) }. R' D9 w% g0 R. t6 n: _
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));3 U- w, A\" Y/ b1 B; g
   
8.    }  ]- q1 Q3 O1 E6 k8 ]6 `5 T8 R\" c+ X
   
9. }

复制代码

Matlab代码：

1. tic0 Z* V: @( U& o- n) \8 Q1 n9 ^
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);$ ?* Y* E* }$ H2 O
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))8 w6 F) \+ x8 T: p
   
4. toc# h: L  L6 V; {+ p: A+ S
   
5. 
   - m8 J7 o; s9 l
6. s =
   2 \% D; d% s7 g, F* Y
7. 
   . y) t+ m4 V* H% o& H: o
8.   1.0086e+0063 A. r3 V8 T+ t+ C5 U6 P7 c. t: [
   
9. 
   0 b# a6 }8 g* B' S. w4 k4 N
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. / ?- o* \: V6 |\\" m9 t  E
4. mvar:
   
5. . N5 h! ]! Q# s% G) D' M\\" W) Y4 C
6. t=clock(),
   
7. 5 t8 d% \4 @$ k; w9 ?+ j' @
8. oo{
   
9. 4 t5 S1 Y2 T& f8 Z. V( b
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 8 g/ A& d6 e2 m\\" m2 m, b: \9 u. 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]
    
13. \\" G+ I4 _* H9 D3 j6 k
14. };5 }! o# ~4 B8 x! l& H9 g
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">/ ?- o* \: V6 |\\" m9 t  E</span> mvar:<br /> <span style=\\"display:none\\">. N5 h! ]! Q# s% G) D' M\\" W) Y4 C</span> t=clock(),<br /> <span style=\\"display:none\\">5 t8 d% \4 @$ k; w9 ?+ j' @</span> oo{<br /> <span style=\\"display:none\\">4 t5 S1 Y2 T& f8 Z. V( b</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">8 g/ A& d6 e2 m\\" m2 m, b: \9 u. U</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\\">\\" G+ I4 _* H9 D3 j6 k</span> };<font class=\\"jammer\\">5 }! o# ~4 B8 x! l& H9 g</font><br /> [clock()-t]/1000;

结果：
: e& h, e, L* e9 M2 C- |1008606.64947441
0.641

- \& b7 S- M3 W$ F0 [Forcal比Matlab稍慢些。

6 R5 Z3 z* F6 Z% x－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic4 e5 a6 ?; L- ~) t, d3 o) ?% ]
   
2. s=0;
   : ]\" V4 N) A, w  K0 {
3. x = 0;4 X% E1 w# a! e4 X\" l5 ]7 V2 \
   
4. y = 1;
   \" b* s6 E$ |% D9 j* r/ Y5 c& Q
5. while x<1   
   9 S8 b7 F  q# m9 a
6.     while y<2;        
   ( {$ q5 j; _5 N! ~$ d* @; R
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
     v+ e& t! i7 V* J6 Q
8.         y = y+0.0009;        
   / E$ h- b8 f  O$ A# [5 Y, R3 c# W3 A
9.     end
   , ^, ]# _0 W1 k- o3 m* h5 W
10.     x = x+0.0011;6 L1 q, I* l7 T0 w! D
    
11.     y = 1;
    6 R6 g1 @9 y: g8 `# W
12. end2 U) B. q- P& ?- r- p2 L1 B
    
13. s( q$ ?\" j- }) R, U/ N
    
14. toc
    : `* G\" i: h: Z
15. \" v1 O5 i$ w8 c1 T) G! Z5 y; m# a) s
    
16. s =
    : K$ S  z: U& V
17. 
    2 Z1 c2 a; u8 {0 i
18.   1.0086e+006
    1 r- m% r5 g* g  Q& z+ l
19. 1 s2 F6 M$ _0 w0 I2 L, z
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:& l0 f8 ^( C0 w
   
2. t=sys::clock();) X% |) z2 _\" h$ |
   
3. s=0,x=0,
     `; a2 d$ k5 v\" X/ r1 i( z
4. while{x<=1, //while循环算法； 5 j5 `% }5 s: }: c3 D2 e7 C# {
   
5.     y=1,
   0 g\" t( k( \( {( Q/ h
6.     while{y<=2,
   . u; E, W- C7 q% Y# v0 O! O: b
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), \" E6 A( E) i0 Z; W9 o2 p
   
8.         y=y+0.0009 ) L: [1 B2 ]3 X. n5 L\" r  N
   
9.     },
   $ }2 ?9 G/ @: ?' n8 i7 U
10.     x=x+0.0011 ' F# i( f$ U% o8 }6 s+ I( [
    
11. }, + n8 o$ \& W; M1 Z4 l4 d  h4 u3 e
    
12. s;, T+ I6 j& V* D0 G5 z
    
13. [sys::clock()-t]/1000;

复制代码

结果：
9 _" N0 |9 v( K- Z) I1008606.64947441
0.734   //时间，秒
6 C$ U4 ?6 X; b& S- {
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
: {+ U) ]: Z4 {; a
& ^" F) R1 H7 _－－－－－－－
- I$ l6 E; z, I* v' V+ F- I6 k
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   . o9 D8 l9 e9 b\" w
2. t=sys::clock();
   7 h9 G4 T& n) W) l8 w\" z; r
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 {% k3 j( ?) C3 U
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   4 B. ~3 q% m5 e* L, v+ N
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
" Q& }1 n6 k% k) A
7 {5 p* K1 H& y* e0 ]9 @/ X- Smatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));, ]0 {  H: P3 }: }1 T% a) J' {& S
   
2. tic
     x4 m0 y1 t) z* V
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);) |3 J+ U& N7 |7 j& c2 ~+ M* L3 Y# b1 o
   
4. sum(sum(arrayfun(f,x,y)))9 k9 l9 t6 G% P% ]; R' u% S5 y
   
5. toc' ]3 f3 I- L; ]\" T: Y8 a$ k
   
6. 1 m$ r! d/ h  ]2 H! n; h
   
7. ans =
   ; d$ t- r) {\" {0 ^& e$ y0 b- v! h9 N
8. : i: z* x8 R( V1 J+ i1 ?
   
9.   1.0086e+006
   8 \\" R5 }) _8 s/ _& P
10. : f- s( [1 a3 y$ d$ P4 p% ~
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. # k2 n4 M8 w; z$ k, W& H& x
4. mvar:
   
5. . p, E! Q( ~9 ^* z7 `) [\\" q
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 8 c8 u4 S; ~8 |7 Y$ {  t
8. t=clock(),
   
9. + H! S\\" H+ E: v3 p4 W/ H  {
10. oo{* X: X/ h( m! c6 j4 W0 v% z
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 0 V  p\\" o, M/ q+ I9 a/ z
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. ! s& W3 b$ h  [
16. };
    
17. + Z- Y6 l7 a/ s8 {$ f& b6 @
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\"># k2 n4 M8 w; z$ k, W& H& x</span> mvar:<br /> <span style=\\"display:none\\">. p, E! Q( ~9 ^* z7 `) [\\" q</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">8 c8 u4 S; ~8 |7 Y$ {  t</span> t=clock(),<br /> <span style=\\"display:none\\">+ H! S\\" H+ E: v3 p4 W/ H  {</span> oo{<font class=\\"jammer\\">* X: X/ h( m! c6 j4 W0 v% z</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">0 V  p\\" o, M/ q+ I9 a/ z</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">! s& W3 b$ h  [</span> };<br /> <span style=\\"display:none\\">+ Z- Y6 l7 a/ s8 {$ f& b6 @</span> [clock()-t]/1000;

结果：
1008606.64947441
0 A% b9 F7 R; h4 b+ m3 H* c- U0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
/ I; g0 P3 S2 V' K) h3 h  Z# z. B, k
－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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