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

forcal

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

1. //用C++代码描述为：) g) F, b4 l4 }& {6 L+ j6 _
   
2. s=0.0;
   ! k1 e2 R& y  h\" \7 _' M4 Z
3. for(x=0.0;x<=1.0;x=x+0.0011)
   2 \  \. G( v/ r  W7 m
4. {
   1 ?! F( w; h0 _8 F% S, W' B
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   ) m1 ?6 W3 d) [
6.    {
   5 T( M7 v6 I6 U
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   7 a. G$ j/ N7 M- G# Z
8.    }
   ) {; M. u* ^0 b4 \0 b9 T9 d4 E8 ^
9. }

复制代码

Matlab代码：

1. tic
   $ V! C8 s2 y( C4 w( A! t
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   . t. J9 k. A0 g2 S0 Y7 x9 {
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ! ]9 j& I, x( }
4. toc
   * _# d6 W# X# \, }, ^* ?
5. 9 F: G3 o\" p2 P+ J) m% q3 u
   
6. s =
   . T3 L5 _: u9 m3 |$ _$ D
7. 2 r; Z6 U9 t\" ~0 {$ e\" S; N# s
   
8.   1.0086e+006( O% [3 j8 A9 P3 }- Q
   
9. 
   0 B, _8 V; p2 q+ X\" @
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];! h  Q\\" G6 Y\\" |
   
4. mvar:9 Z* t0 B: }0 l( E6 y
   
6. t=clock(),
   
7. $ t2 h) ^# B) j5 P  l
8. oo{- k\\" O& N\\" R6 T1 X: p
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. - [% Z# F6 k- }0 z9 H& V& J
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]$ T* f9 Z& ]9 j3 k
    
14. };
    
15. ) l5 I. n4 }& ]9 B# l
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">! h  Q\\" G6 Y\\" |</font><br /> mvar:<font class=\\"jammer\\">9 Z* t0 B: }0 l( E6 y</font><br /> t=clock(),<br /> <span style=\\"display:none\\">$ t2 h) ^# B) j5 P  l</span> oo{<font class=\\"jammer\\">- k\\" O& N\\" R6 T1 X: p</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">- [% Z# F6 k- }0 z9 H& V& J</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]<font class=\\"jammer\\">$ T* f9 Z& ]9 j3 k</font><br /> };<br /> <span style=\\"display:none\\">) l5 I. n4 }& ]9 B# l</span> [clock()-t]/1000;

结果：
! Z/ m! m( X& b9 x1008606.64947441
0.641
# B, D/ V) a, p) @9 a
7 G$ G5 h0 A/ `' G) r# d( EForcal比Matlab稍慢些。

－－－－－－－－－－
; d/ m1 g. a% k
再看循环效率。
1 M! D) z6 O1 m0 u/ S2 J5 |
Matlab代码：

1. tic
   2 c& C! q' K+ E8 G% W7 O& F
2. s=0;$ @3 n$ P0 u6 _; h  s4 k
   
3. x = 0;6 h* L# w9 D  t0 e
   
4. y = 1;
   1 x( b- L7 c: z. U- m- L
5. while x<1    , P2 q0 q- U3 K8 U) l' u8 U: X
   
6.     while y<2;        
   5 x% ^6 ~+ Y5 X3 |- D
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 i\" a0 Y& x) {
8.         y = y+0.0009;        
   ; H\" v) K2 d8 y5 I: [5 R9 K
9.     end
   1 s$ u+ c! w! C3 s* L! h; d* I
10.     x = x+0.0011;& w& H' T) G. l+ F0 |* V/ v, E
    
11.     y = 1;
    / P1 e% `\" u9 v
12. end
    ' r4 Z% Z\" p2 n/ g) i
13. s
    0 `- o# v0 I; @; d. L
14. toc. L: O9 m' f$ B7 w) M( R1 d# E
    
15. ( S2 Q8 w# _5 p( C% v1 J1 V2 n
    
16. s =
    / o; \) j; j2 h% l
17. + l4 t  q# g4 R- {' u1 b, p  Z
    
18.   1.0086e+006
    ' _% _/ b- Q9 ?9 M4 s. ~4 g* D
19. 6 t7 ?( \( N8 p4 Z- A/ D2 G4 p) M
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   ( Z( b% M$ E; b
2. t=sys::clock();
   $ W4 X; u3 U# ~, s+ m2 O
3. s=0,x=0, 6 f% l# d$ U( T3 q5 m
   
4. while{x<=1, //while循环算法； - d# S; q3 E4 |$ J3 \
   
5.     y=1,
   ( }' c. U5 s2 I1 t. ^, H3 @) }
6.     while{y<=2, & [- @& J/ C+ k4 {, x# P- ~
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ; F# `% {! v% g* X4 g' z( `
8.         y=y+0.0009
   9 c$ m+ o; N+ L% ?3 Q2 ~
9.     }, + ]+ a: w+ S8 V' Q- u& x1 K
   
10.     x=x+0.0011 % l/ ?& h/ {9 E# w2 N4 V
    
11. }, 1 A7 F/ @& @# c# L( o4 K
    
12. s;
    , Q/ a6 }3 @+ J$ V3 q/ W
13. [sys::clock()-t]/1000;

复制代码

结果：
8 N7 |5 Z- c' ^( u: q0 o: b1008606.64947441
- J+ b% u/ ]7 N7 q, k( R0.734   //时间，秒
: B* Y. b3 @6 Y. _( L0 G. g
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

- f+ ~1 P# Q" X+ r－－－－－－－

/ A+ E. q. z8 j* wForcal中还有一个函数sum专门进行这种计算：

1. mvar:* ~2 i& \\" U& c6 u\" l; @* A
   
2. t=sys::clock();2 j- |. m; M. @
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 7 j5 V6 y1 X8 w) O8 R9 W- r1 B
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   9 B1 g' [% P& l% X) A
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
) F1 P+ q# y5 q
: Y& ~7 h- v9 @# @' a1 kmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));7 {& M9 d% l* d2 L3 h\" {! i) j\" H) x
   
2. tic# h) W1 h0 T% F- N7 T
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ) y2 B3 K, h* N\" U$ X( U
4. sum(sum(arrayfun(f,x,y)))
   \" h/ S+ _! X/ ~0 \. N; M, M
5. toc
   9 v! l( M\" P  Y
6. 
   + w* z. c$ ~% f$ C7 K
7. ans =
   ! n  O+ }5 v; n1 Q
8. 4 a* |' Z! y. E1 k. z
   
9.   1.0086e+006
   8 f2 b\" F. O\" _% T+ s
10. 
    + K. y( q; G9 A\" x$ z
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];6 |2 I! h1 k. m9 K7 @1 I# a( L
   
4. mvar:# Y( `/ O) {: g/ {6 q
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); / f# Q$ N2 l: S) M$ `
   
8. t=clock(),3 c  [# c  K& A! L4 v
   
10. oo{
    
11. # z0 F& Q% o6 Y
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 5 I% Y  P# a0 I5 Y% @, n: e
14.     arrayfun[HFor("f"),x,y].Sum[]3 z- u4 J1 b6 T, k8 b2 L
    
16. };
    
17. : {: z7 M3 h- {1 D
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">6 |2 I! h1 k. m9 K7 @1 I# a( L</font><br /> mvar:<font class=\\"jammer\\"># Y( `/ O) {: g/ {6 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\\">/ f# Q$ N2 l: S) M$ `</font><br /> t=clock(),<font class=\\"jammer\\">3 c  [# c  K& A! L4 v</font><br /> oo{<br /> <span style=\\"display:none\\"># z0 F& Q% o6 Y</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">5 I% Y  P# a0 I5 Y% @, n: e</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">3 z- u4 J1 b6 T, k8 b2 L</font><br /> };<br /> <span style=\\"display:none\\">: {: z7 M3 h- {1 D</span> [clock()-t]/1000;

结果：
1008606.64947441
% G, d* S' s( J+ R6 e1 A: Y# c2 e0.735  秒

; m# a; s# S* m9 j2 m/ r可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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