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

forcal

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

1. //用C++代码描述为：
   , J/ m2 g) l$ Y. o
2. s=0.0;
   / w\" z0 O$ w! ~7 p; n8 y
3. for(x=0.0;x<=1.0;x=x+0.0011) . f- H( i\" W! e9 q4 z- h/ ]$ q, j. b
   
4. {
   . B  Q1 z8 f9 S8 K
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   4 I% ~1 m2 q% k9 V
6.    {
   5 N) @- g, M, G. I: _4 U# _7 b! f
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! e: }\" ?' f2 e# w. r1 w7 X- G
8.    }& l( t+ h\" Z& J7 e. @9 c
   
9. }

复制代码

Matlab代码：

1. tic
   ( F* B5 m, o0 N, Y  a
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);  m7 w2 l- b% C$ ]6 A7 y
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   2 s6 r3 V- i0 s) B1 Z
4. toc0 A2 n0 {- `5 {( n6 [
   
5. 
   + P, ~3 P: P* O3 H
6. s =! x  M+ m6 m7 c( `6 e: ^, }* P' W
   
7. 
   & T& x9 M. ?! w. k8 t8 j* }3 h1 c
8.   1.0086e+006
   / b\" u  ^0 D* U6 X% T, w\" @* ]0 ^. s- x
9. ( I3 Z  S% e, `* D  L& r
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 5 l9 N8 P3 @7 O# U+ q' n9 T
4. mvar:
   
5. ' h1 e1 D* [: A1 u# u/ D- B: o
6. t=clock(),
   
7. # g* M. s5 x' b: e4 q
8. oo{* N: `$ [4 a9 ]* w/ l6 |  C; }4 V+ V
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],; M$ y$ ^! R+ o6 Y& Y: B/ w, \
    
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. 9 B% o5 v/ A) c  l7 L* j
14. };9 ^# w) f- u8 O# Z
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">5 l9 N8 P3 @7 O# U+ q' n9 T</span> mvar:<br /> <span style=\\"display:none\\">' h1 e1 D* [: A1 u# u/ D- B: o</span> t=clock(),<br /> <span style=\\"display:none\\"># g* M. s5 x' b: e4 q</span> oo{<font class=\\"jammer\\">* N: `$ [4 a9 ]* w/ l6 |  C; }4 V+ V</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">; M$ y$ ^! R+ o6 Y& Y: B/ w, \</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\\">9 B% o5 v/ A) c  l7 L* j</span> };<font class=\\"jammer\\">9 ^# w) f- u8 O# Z</font><br /> [clock()-t]/1000;

结果：
! x* U7 [! H/ Y, o  l) i1008606.64947441
# f3 h6 f5 a% m: n4 _. q9 S7 J) U0.641
: ]1 h9 l$ H0 p0 i# Y, S6 l+ ~4 S
) ~4 ]/ [' n% d9 ZForcal比Matlab稍慢些。

. v3 @7 E: s7 k! R+ T. S8 r－－－－－－－－－－
) t7 f- [4 G1 V, ^+ ^  K7 k
再看循环效率。
% _* V) F3 A* U& {9 h: ]# y# a
Matlab代码：

1. tic( D3 \& i1 e( h
   
2. s=0;
   5 R: e' W; K4 Y
3. x = 0;
   - g+ q/ g) v  W/ U. }
4. y = 1;% B' K& W; }5 `9 w
   
5. while x<1    3 ]5 J/ v1 E2 J* V/ l/ w1 A
   
6.     while y<2;        4 c% M7 d; r+ Q\" a
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));\" @1 s% |/ Q* p- G
   
8.         y = y+0.0009;        
   $ X; S! `& \$ V$ V
9.     end
   - p  \) R1 d0 H6 e1 T1 C1 @2 N
10.     x = x+0.0011;; P  _* g* I4 O9 v
    
11.     y = 1;: w\" s& u- I5 Q: W\" c& g
    
12. end
    ( q7 {7 x' i) o\" ~
13. s6 l' ]1 u5 d& J\" ]; x3 X; E
    
14. toc  R5 }, J& G  v; V& h9 |. o9 I
    
15. 
    5 L% |/ M( g! H- T* c9 k
16. s =
    7 `- s* t: C' |  ]
17. + j3 P8 v$ m* R  p
    
18.   1.0086e+0065 E2 \$ o8 R9 C\" K' o2 d
    
19. \" n/ b7 F. f; A
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:2 R0 g! E; p2 O, r! y- r
   
2. t=sys::clock();
   * _6 M- p/ u. y! J% ~
3. s=0,x=0, 8 r! R( C\" O: O3 N3 d/ Z
   
4. while{x<=1, //while循环算法； , r: g) s* A, D\" Y) V  @+ T
   
5.     y=1, - g$ W3 a- Q0 j: r/ A4 M7 D
   
6.     while{y<=2, 0 C' a5 w8 l2 b+ ?1 S/ J; a, r3 ?
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 5 v+ Q  k7 j7 x# C  H: [1 o
   
8.         y=y+0.0009 & H; E' w$ Y/ [9 z
   
9.     },
   5 S9 Y3 J/ R9 n8 U+ d  J+ }% d9 I- m
10.     x=x+0.0011
    6 v3 y- S% e5 J$ r
11. },
    ( R3 h+ M0 p* ~  `
12. s;
    1 r5 s9 \' N1 m  D3 @+ j
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

  H$ T' ~1 V# S( K$ K4 m: E我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

2 W: G1 B' T: G, q* |－－－－－－－

  O4 ^  n. B6 v6 n: CForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   ) o3 J- ^# L) j+ [6 T
2. t=sys::clock();, U; t& G' Z2 ]$ }, f% K
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 0 w( p! Z: y) A\" S8 X- C- d1 }
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   5 W( @8 U3 X/ z- s2 ]
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
5 @5 L$ g+ y: P9 v0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

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

+ E4 T* j! g0 ]4 C- d6 omatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 n7 ?/ \  k' H
2. tic; b. I) c5 A+ l* s+ F1 ?
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   $ H$ g9 R' Q! }
4. sum(sum(arrayfun(f,x,y)))( t9 _! `5 X8 z+ P# u
   
5. toc$ N& l! s2 T; z- Q' T3 R1 q/ U
   
6. 
   8 }/ C# p! d; x# x; x
7. ans =! y' F, D( k( ^' y) [
   
8. 
   1 k/ Q  d& r) b' c4 {
9.   1.0086e+006
   + V$ _\" }* E4 l4 e* G
10. 
    , v6 ]7 e% E. V! j+ T0 `\" C! |+ D
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. - V\\" A\\" z( T+ \/ d/ ]/ i
4. mvar:\\" C) Q1 }6 j5 J0 K; Y\\" V1 q
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 5 L1 ?7 D* t- ^! F' S
   
8. t=clock(),3 h$ t# D! J9 c. C0 e! V
   
10. oo{) F4 }- f. j) o! o' G4 ?3 c
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],# q\\" A\\" M9 {/ I! A
    
14.     arrayfun[HFor("f"),x,y].Sum[]& Y& J. p; U5 L  A4 A# x( z
    
16. };3 M  |1 S3 ?# c. a9 J0 C
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">- V\\" A\\" z( T+ \/ d/ ]/ i</span> mvar:<font class=\\"jammer\\">\\" C) Q1 }6 j5 J0 K; Y\\" V1 q</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">5 L1 ?7 D* t- ^! F' S</font><br /> t=clock(),<font class=\\"jammer\\">3 h$ t# D! J9 c. C0 e! V</font><br /> oo{<font class=\\"jammer\\">) F4 }- f. j) o! o' G4 ?3 c</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\"># q\\" A\\" M9 {/ I! A</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">& Y& J. p; U5 L  A4 A# x( z</font><br /> };<font class=\\"jammer\\">3 M  |1 S3 ?# c. a9 J0 C</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
8 O; ~/ z4 @$ L
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

5 _* J3 x! F+ \) F% Q" x从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。