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

forcal

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

1. //用C++代码描述为：# F. J2 c# O2 Q2 F3 I
   
2. s=0.0;
   $ P# A2 f& k% x( w: m: D
3. for(x=0.0;x<=1.0;x=x+0.0011)
   . R: c6 ~& J  B; q# S
4. {
   0 S8 N# k5 j0 o! P6 J1 q4 R) ?
5.    for(y=1.0;y<=2.0;y=y+0.0009)' n8 s1 v5 C! L& @% a$ M; g; E
   
6.    {5 q* c  \7 u% S: _6 ~: ~
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ v# A/ Y  P# P1 s/ s
   
8.    }
   ; r/ H3 A+ k9 G9 }9 f3 ?
9. }

复制代码

Matlab代码：

1. tic) x2 B7 H, Y6 P5 L5 C* H
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   8 J% v( `, K: f# m5 |
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))\" `/ }$ N; q% A( H
   
4. toc  O! z* L3 V' v* d- P5 k0 n
   
5. 
   % w, y; r+ i8 j! |) i# \
6. s =8 `# c) e; O0 A' i* s* R6 ]% T  k
   
7. 
   5 x\" h! ]# o7 A! V% w( k) D/ s7 `
8.   1.0086e+006
   : \\" z, K' X; S9 J
9. 
   6 Q( {$ W\" r1 B6 A5 g( b. a
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];( g' s) N5 h5 i/ L  l$ T% C. f
   
4. mvar:
   
5. - U- a& O! E1 {6 d0 A
6. t=clock(),1 N# g' U9 K8 h/ _\\" m4 {% b9 U
   
8. oo{( h; b3 z  q  O3 v2 s
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. ) q2 ?. @0 I8 F
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. ) j2 ~3 F\\" G% \. G7 Q' ^; f
14. };
    
15. $ f: y# `  O6 Y7 X( v
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">( g' s) N5 h5 i/ L  l$ T% C. f</font><br /> mvar:<br /> <span style=\\"display:none\\">- U- a& O! E1 {6 d0 A</span> t=clock(),<font class=\\"jammer\\">1 N# g' U9 K8 h/ _\\" m4 {% b9 U</font><br /> oo{<font class=\\"jammer\\">( h; b3 z  q  O3 v2 s</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">) q2 ?. @0 I8 F</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\\">) j2 ~3 F\\" G% \. G7 Q' ^; f</span> };<br /> <span style=\\"display:none\\">$ f: y# `  O6 Y7 X( v</span> [clock()-t]/1000;

结果：
1 N3 z. o0 s2 {1008606.64947441
0.641

4 A& a, f$ E  m8 VForcal比Matlab稍慢些。
$ [; S. G" u* R% U& @: x3 p
0 @2 n* e7 V7 Y# y－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic6 u\" j0 j1 g7 R0 m6 S
   
2. s=0;% }\" d- `: ?1 J* q: e0 @! |
   
3. x = 0;
   2 E+ ~4 u+ ]# M3 n. R  K
4. y = 1;% ?9 `1 u* o# g8 i
   
5. while x<1   
   & w, |3 e+ [( f: @. }% h. j' I
6.     while y<2;        
   % Q$ A; k; V, ~8 n( J: T  t
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - _8 U& C\" \5 W  J( g1 Y+ r
8.         y = y+0.0009;        / G1 F/ b8 H$ g5 _: V# z
   
9.     end
   , J+ v4 d3 Z/ m0 b8 o0 ~/ V
10.     x = x+0.0011;! T5 K/ w* ^0 R\" @
    
11.     y = 1;+ x6 i+ u/ _& T; K/ ]. K
    
12. end
    - n& @4 P  j: l  _
13. s\" ^% F. b0 f4 H, j2 }# |, F. g
    
14. toc
    7 P, C4 x8 o! S
15. ( W4 U; l  e5 b/ O- o1 [0 O
    
16. s =! P3 H1 X2 N4 ?7 l* {: Y5 p
    
17. $ T  ~5 `3 W9 Q\" T- ?
    
18.   1.0086e+006
    : M: m6 y2 w5 _% }, q7 ]
19. 
    \" e9 n$ Q6 B+ _\" O- m0 G* p( l
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:6 t6 v# ~& e% a
   
2. t=sys::clock();# X  z; {3 S1 a$ F+ x  q: Z
   
3. s=0,x=0,
   , w\" M' H9 H$ n* ?0 }9 S
4. while{x<=1, //while循环算法；
   ! T' s0 K0 b; F( q1 Q. h# m1 ^( i
5.     y=1, 9 q& ~6 g* Z3 r: P5 u
   
6.     while{y<=2, ( E0 S, I, u& \& @# T( Q9 V2 j
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), - ]4 m8 A' ]0 g3 u: s
   
8.         y=y+0.0009
   2 c* v& S% B3 K' `' ^$ p
9.     }, 1 I- q\" W. Y' ]/ I: }8 o; h
   
10.     x=x+0.0011 / e+ q( [\" _2 X; ~# R
    
11. }, : B3 W  U2 Y6 L- y9 A3 t
    
12. s;
    . v$ D- S$ Q. r\" [5 m% M0 j2 [/ U
13. [sys::clock()-t]/1000;

复制代码

结果：
6 R( q- }+ ~5 u6 k; ?1008606.64947441
: a$ F* e8 D5 c3 v) A$ T8 E! t0.734   //时间，秒
. Z5 d- v+ Y9 d" Y
  Y, s' E$ y3 ^; i# L我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

; G1 A1 A& l' B+ h# V5 C' n－－－－－－－
) g9 E: G$ `$ Y7 o
* r! R6 n9 k" c7 h% S3 {3 OForcal中还有一个函数sum专门进行这种计算：

1. mvar:. D3 n6 I2 d% F1 t6 c  ?
   
2. t=sys::clock();
   . R/ {4 \1 Y\" `/ F$ l4 B. N, \* B
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 9 V3 f: y. X\" Y# A
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];6 p7 V3 W9 Y3 B7 I4 C5 k
   
5. [sys::clock()-t]/1000;

复制代码

结果：
7 h0 f! r, N1 \: @, X; v$ U2 L! c1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   6 W8 D4 L  T% i+ u& g( ~6 ^7 V( [
2. tic
   7 _- {7 ^: J7 m8 ?+ B/ Q
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   \" r4 O3 D) k7 y% k/ P% [# A
4. sum(sum(arrayfun(f,x,y)))- j& _9 t$ y: s6 I* ]+ j
   
5. toc7 J+ E# l+ N6 Q6 s4 O9 L# v3 j
   
6. 
   + Z6 c  ^# t3 U3 Q$ H4 k3 U# K
7. ans =: y7 Y7 Q5 z% g
   
8. ' |9 c* H# }5 K\" m/ R
   
9.   1.0086e+006$ M\" ]% E; R7 z: l' {- e
   
10. 
    ; d6 q9 {( F  @5 G% R' ]
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 6 y$ [, Y/ x% e. c8 H
4. mvar:
   
5. ; W9 N; l* v4 p
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. . G% z9 G  M: J( g# A6 `
8. t=clock(),
   
9. , R8 d/ v# \6 h+ ^1 c3 u9 }
10. oo{( g4 K' H/ C6 q5 R5 j
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. $ k' O$ x' t! U  f
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 0 U  r9 i/ v! N* f
16. };
    
17. 7 P# i& H; Q2 \2 |- o1 L4 y
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">6 y$ [, Y/ x% e. c8 H</span> mvar:<br /> <span style=\\"display:none\\">; W9 N; l* v4 p</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\\">. G% z9 G  M: J( g# A6 `</span> t=clock(),<br /> <span style=\\"display:none\\">, R8 d/ v# \6 h+ ^1 c3 u9 }</span> oo{<font class=\\"jammer\\">( g4 K' H/ C6 q5 R5 j</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">$ k' O$ x' t! U  f</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">0 U  r9 i/ v! N* f</span> };<br /> <span style=\\"display:none\\">7 P# i& H; Q2 \2 |- o1 L4 y</span> [clock()-t]/1000;

结果：
: ?1 p9 E* b$ R0 I* ]! R1008606.64947441
1 b: E0 a2 W; ?# W+ m# j# U- K4 w0.735  秒
7 a' _- l/ I# S' h
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
# \' P4 F" _$ n9 C
0 y; `5 |% d- O+ c: E9 R% F－－－－－－－－

+ P0 I6 d* v+ |: O4 b" b/ p1 b* c' T从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。