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

forcal

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

1. //用C++代码描述为：9 X, g8 k' s6 I. D/ q1 U# ~: ~% `. t
   
2. s=0.0;
   ! W7 m* f  o& D4 z5 g) S* w
3. for(x=0.0;x<=1.0;x=x+0.0011)
   / B4 a1 S9 t8 \+ [( O: H
4. {\" V7 w7 p* P- `1 ]) i: q( F
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)* ~5 e\" ~# U* {3 `0 [7 N6 @& K, q
   
6.    {. v+ g1 y! V7 i7 Y! z3 Q
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));8 K3 r4 U+ x5 H& Y; [
   
8.    }* |7 ?2 V& ~$ n. l, B5 z
   
9. }

复制代码

Matlab代码：

1. tic$ v\" U% i1 V$ `- L) N2 c
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);0 h: \! R% l# c+ e5 U) l( I2 s* e
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   * @  O* \. q; v& B- B
4. toc
   / S4 L\" H; p' E; e/ l* v
5. 
   ) N$ W4 }9 k4 a$ f) y\" R
6. s =
   ' D- ?5 |\" \& w  F/ v
7. ; l1 V4 t\" |1 z/ l2 C4 q& {
   
8.   1.0086e+006
   / Z! ^7 |* L$ h1 F$ |# M# b: d4 ~/ l' v
9. ! P  a' R4 E3 P
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 5 Z  o$ _7 h3 N2 t; n0 F; ]$ T& z; X
4. mvar:4 [5 i' [! p1 s$ R: M: D6 Z
   
6. t=clock(),
   
7. # s& T: E$ K( _$ G8 C: l: U$ M
8. oo{
   
9. 7 k\\" q7 c7 B3 J
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],# e8 M- D: i$ P3 r5 ?' O
    
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. 4 Z; S/ S8 c* B' X' i
14. };: G% }* c3 z  t: O' M& N' ]
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">5 Z  o$ _7 h3 N2 t; n0 F; ]$ T& z; X</span> mvar:<font class=\\"jammer\\">4 [5 i' [! p1 s$ R: M: D6 Z</font><br /> t=clock(),<br /> <span style=\\"display:none\\"># s& T: E$ K( _$ G8 C: l: U$ M</span> oo{<br /> <span style=\\"display:none\\">7 k\\" q7 c7 B3 J</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\"># e8 M- D: i$ P3 r5 ?' O</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\\">4 Z; S/ S8 c* B' X' i</span> };<font class=\\"jammer\\">: G% }* c3 z  t: O' M& N' ]</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
6 l- r# z0 K& a' F4 n) F0.641
7 L9 _5 v, [5 ]6 G/ i& x! P9 x1 n
Forcal比Matlab稍慢些。

& W  p! [0 s, g: e' W$ v－－－－－－－－－－
% {) H2 I# G% a
再看循环效率。

" L7 T. R  V1 i0 B3 y4 C0 xMatlab代码：

1. tic
   . @1 Q/ R# h' A5 K0 J; k# ~
2. s=0;- K+ @- P1 n  U$ L5 A
   
3. x = 0;
   % Z; m8 Z4 `2 q; F( G
4. y = 1;
   . P' [1 ?! _/ B2 E0 c
5. while x<1   
   7 a! v# A$ i! [  h- d6 C) O
6.     while y<2;        2 e$ z4 Q: C/ p
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % v1 y5 k! t; T* E' N
8.         y = y+0.0009;        
   ! o* M  b% }8 i7 j: N
9.     end
   7 h0 a) _  l5 l- c6 }
10.     x = x+0.0011;
    3 B: _% l5 ~1 L  H/ q3 A6 y9 `
11.     y = 1;3 O- E0 X$ ^* N
    
12. end
    ' J8 E' k3 z, L. n
13. s* a\" \0 u! @: _6 e) e
    
14. toc( s  X& |0 M, z$ b) e
    
15. 0 f! c* s; f0 w. j/ J! `' f
    
16. s =! j; b0 G- Q! U; K
    
17. 
    * Z# X, Y) V; a6 b' }$ l
18.   1.0086e+006
    ) H* m+ c8 p! e0 |# A* n
19. 9 w4 k/ J3 v9 m3 E0 ~
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   7 l( y* e: D* f% N; [, m5 r
2. t=sys::clock();
   $ w/ K/ m7 }3 a5 Q
3. s=0,x=0,
   % v1 O5 z8 a' O1 o
4. while{x<=1, //while循环算法；
     f3 P' P\" |( W
5.     y=1, 1 h\" |9 J& @8 G; w' c; k\" J
   
6.     while{y<=2,
   . p. q1 @- N* D% b\" X+ T% ]/ g- m
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   6 O3 w' b\" Q2 g) W5 r
8.         y=y+0.0009 ; i% k% Q: s4 Q8 Z6 t
   
9.     }, 8 w\" m* V0 W* s! p! C' E8 Z0 z. X
   
10.     x=x+0.0011
    9 m4 d3 ?% F) v9 t) P  e
11. }, ( M2 t7 _7 h4 y! @5 c  m
    
12. s;) A% w( F- |: `0 g2 R$ @
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
* r# s5 h0 `! q2 Z0.734   //时间，秒

8 _1 W9 w, V$ @我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   6 H3 p. @6 N9 B* m$ v, n( Q  U
2. t=sys::clock();4 g5 d/ n9 v, m1 ]$ G  R9 T
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ! V8 e7 ^$ b# z\" h
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   - l\" ^( D+ h+ O4 @
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
4 Y( H4 R' p, W$ ?1 _7 _  x( w0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
; c6 r- I0 N2 n2 F/ j1 E
* [* }* j4 p1 H) K0 Qmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));7 u  h( @' O: t0 u. U3 k  |
   
2. tic4 n8 k; @: Z' `/ J, G# Y6 x2 D# n$ i
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);* C# V) R( a& J7 e5 N
   
4. sum(sum(arrayfun(f,x,y)))' K/ t0 j: d5 }$ `* g+ q/ R
   
5. toc: X' O* {0 I* b& f\" ?6 \
   
6. 
   ! |9 u1 s\" ?1 w4 p5 F
7. ans =
     u/ c8 m% E4 Z\" m7 o& h( \
8. 
   ! e: T- V6 l. }8 Y
9.   1.0086e+006- }0 Q2 l: Y- m
   
10. 
    7 g4 r; w\" g2 D
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];0 B2 T$ g\\" E. V6 D
   
4. mvar:
   
5. * F& R  w8 ?9 o  k* ]% ]. j- Y. i0 \- K
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. # U; O0 m& W8 n' {\\" f6 Z
8. t=clock(),
   
9. \\" t7 O% y7 N* r/ @) V
10. oo{9 S, I' c# B# h5 x# q' S
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 1 f/ _\\" }, `% D2 f- L. O0 z
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 5 L2 ]( i( S. ?5 S: f
16. };  c2 M+ O9 Y7 W$ x8 X0 A, V\\" ~
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">0 B2 T$ g\\" E. V6 D</font><br /> mvar:<br /> <span style=\\"display:none\\">* F& R  w8 ?9 o  k* ]% ]. j- Y. i0 \- K</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\\"># U; O0 m& W8 n' {\\" f6 Z</span> t=clock(),<br /> <span style=\\"display:none\\">\\" t7 O% y7 N* r/ @) V</span> oo{<font class=\\"jammer\\">9 S, I' c# B# h5 x# q' S</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">1 f/ _\\" }, `% D2 f- L. O0 z</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">5 L2 ]( i( S. ?5 S: f</span> };<font class=\\"jammer\\">  c2 M+ O9 Y7 W$ x8 X0 A, V\\" ~</font><br /> [clock()-t]/1000;

结果：
9 e( ^1 y1 S9 Z0 `2 B: T/ d! U1008606.64947441
0.735  秒
% L( D! B- v( R2 c; F( W+ o
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

0 L$ _0 ]3 F" k1 \- ~) c9 l8 b－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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