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

forcal

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

1. //用C++代码描述为：
   7 F1 ]$ E3 e6 d: ~, _6 Y5 q' n0 G
2. s=0.0;
   + G, v5 ^* n4 H% \# B4 c
3. for(x=0.0;x<=1.0;x=x+0.0011) + T$ a3 {\" s3 `\" [
   
4. {; U0 L0 x0 W3 R0 ^0 T; `1 `  R* I
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   7 o- P8 F; q\" u/ ~5 I  q! i  z
6.    {
   3 d) o# B- P+ H8 M& W3 k$ w
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 ]3 f  k* s& G6 U
8.    }6 T+ l2 H+ q* x, K6 e
   
9. }

复制代码

Matlab代码：

1. tic6 v* D+ J6 Z' q! u5 o
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);) n/ p6 ?; N# Q3 p9 m* n
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   9 S* Y0 X  N4 N& b* ?- \
4. toc$ R/ }2 _5 ]) X# Y+ u% K: c
   
5. 
   . N  q9 ?3 Q6 V# |* |
6. s =
   8 Z6 h2 \8 n$ a& P* x4 n
7. # d4 g: P2 I8 r: K+ T: M! L8 q
   
8.   1.0086e+006
   ' f4 u) [3 Q; B1 A' N/ d
9. 
   ' G3 {\" H\" t0 E; p, Q\" Z+ D, V: L4 n
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ( O& W3 `. p! [2 ~
4. mvar:+ [( Y# a! n% i; m, |% g
   
6. t=clock(),
   
7. \\" h8 V# L( H9 k. ^. r! R1 ~9 u4 }
8. oo{
   
9. / v  _9 q# a2 \  Z
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: W8 n6 L  v5 M
    
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. ; E2 _8 _4 ~7 b8 r- G
14. };3 R+ [6 V4 ~3 J0 I9 T8 J% Z
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">( O& W3 `. p! [2 ~</span> mvar:<font class=\\"jammer\\">+ [( Y# a! n% i; m, |% g</font><br /> t=clock(),<br /> <span style=\\"display:none\\">\\" h8 V# L( H9 k. ^. r! R1 ~9 u4 }</span> oo{<br /> <span style=\\"display:none\\">/ v  _9 q# a2 \  Z</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: W8 n6 L  v5 M</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\\">; E2 _8 _4 ~7 b8 r- G</span> };<font class=\\"jammer\\">3 R+ [6 V4 ~3 J0 I9 T8 J% Z</font><br /> [clock()-t]/1000;

结果：
: K" |) ]6 j" ~4 s4 o6 f% H1008606.64947441
0.641
+ H; c! `( f6 ^9 S8 {* q( F8 Q
Forcal比Matlab稍慢些。

－－－－－－－－－－
& P2 ?, g# s) S& Z
再看循环效率。
/ r4 X% K  j! {( ?' i0 w
Matlab代码：

1. tic2 Q6 s! e$ G\" n/ H7 T
   
2. s=0;
   8 J8 ^2 P# s. |6 V$ `
3. x = 0;
   0 q5 d* ]0 H, g: H
4. y = 1;
   8 a  c) T, G: `$ c+ u8 n) \
5. while x<1   
   ' u* W4 i6 `; a\" @- q
6.     while y<2;        
   ) c$ q1 M0 \0 q) f  y5 u. X
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 l+ v) ?  Z6 O\" y* Z1 P
8.         y = y+0.0009;        & |6 R: _5 K; t4 Y( ?\" q
   
9.     end
   . ?% G% u1 {$ b$ t2 |9 `% p; u
10.     x = x+0.0011;& Y+ R4 ]1 W) X3 M; s1 B/ V* z
    
11.     y = 1;
    9 z1 A* N* ~) n# ~+ i
12. end+ l- Z2 `, m  P4 e1 M+ ~. y( ]8 S
    
13. s6 M6 P- m' W* L0 d$ b4 p5 R: m
    
14. toc\" X% _+ u\" H* `6 Q* B
    
15. ; F' N5 v6 \* t$ g, v( q
    
16. s =
    8 t  c$ A* ^3 w8 d+ W
17. 
    - |\" |% z\" k' ^  n* e! g
18.   1.0086e+006
    3 {) ]5 Z2 o8 @
19. 
      |' z* O+ O& e
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   & s* `& r9 D' r1 Z
2. t=sys::clock();! X/ O  J8 @4 |\" w
   
3. s=0,x=0,
   6 Q8 d  z' S/ i% ]& U
4. while{x<=1, //while循环算法；
   * q& c& n( f: b
5.     y=1, \" X: h, q% L5 g' W- _0 A6 }0 r  S8 v
   
6.     while{y<=2,
   7 i! c7 K; y! |6 h+ \, L, x
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 0 b8 R& ?' G& j' {& p
   
8.         y=y+0.0009
   $ V# V5 T' p) }3 O  L, w3 s
9.     }, 8 C& ^; g+ Q, x8 e7 W/ U: m1 C# L
   
10.     x=x+0.0011
    ' J( O5 i4 G$ h5 T; @$ H) x9 Q: M
11. },
    9 E- i- ~7 w+ t' M2 d* f$ R. G7 X6 y
12. s;0 S4 \\" h/ }% k% Y: }7 ^
    
13. [sys::clock()-t]/1000;

复制代码

结果：
) |5 t, J, X6 H6 @$ h0 I9 g1008606.64947441
0.734   //时间，秒

1 S; G- Q' y  G0 i, s; l# @/ Z! }我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
, E+ b3 l) Q$ `  Q9 C3 x/ i: N6 e+ O
- R7 C7 q' f( G! N- CForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   2 T+ _5 x& [' N% l  U- l0 Z
2. t=sys::clock();
   ' B2 s5 g1 J, P- {4 H: ]. S) _8 u8 u
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( ]* m3 N\" n0 x! M5 M# ]* J$ W5 d) P. @
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   7 E- z\" w/ p0 g$ ~( ]& J+ A
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

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

0 R5 u+ ]+ T, ?# Kmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));. y$ }9 r% e# e5 M! l  ]# ]
   
2. tic
   + ?+ k. C( m2 P0 @. G
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * `% d) N\" X# A4 C! u& A# u% `
4. sum(sum(arrayfun(f,x,y)))
   \" O* O; X! G# W- P! }0 y! q
5. toc9 H9 Q& I) u4 v- M. X' i
   
6. 3 d  u! W( x7 c7 j
   
7. ans =. `8 A5 b/ D8 k5 \7 R
   
8. \" u6 F. z) p4 j8 }7 Q3 G( ^' R/ r6 d\" P
   
9.   1.0086e+0065 t( U  H! u/ b) T% y0 F
   
10. 
    8 K/ g- A8 W: F- ^- d
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];: k. W\\" B  \0 K$ b\\" n! f4 B7 m* \
   
4. mvar:
   
5. 1 r) S( |) [7 S6 ?
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); / k5 J9 x3 w+ J# c$ J
   
8. t=clock(),* m9 I, ]; \$ n& e* ]\\" Q
   
10. oo{/ I  w6 m8 B: I& }
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],) F/ [' r3 h9 G
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. : G2 l0 Z/ D* F) i+ Z* i( E
16. };
    
17. 8 {\\" |( y' m1 U9 s7 R/ g/ ?; l$ R
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">: k. W\\" B  \0 K$ b\\" n! f4 B7 m* \</font><br /> mvar:<br /> <span style=\\"display:none\\">1 r) S( |) [7 S6 ?</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">/ k5 J9 x3 w+ J# c$ J</font><br /> t=clock(),<font class=\\"jammer\\">* m9 I, ]; \$ n& e* ]\\" Q</font><br /> oo{<font class=\\"jammer\\">/ I  w6 m8 B: I& }</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">) F/ [' r3 h9 G</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">: G2 l0 Z/ D* F) i+ Z* i( E</span> };<br /> <span style=\\"display:none\\">8 {\\" |( y' m1 U9 s7 R/ g/ ?; l$ R</span> [clock()-t]/1000;

结果：
$ k; d7 M0 X3 p" l/ U+ n5 b1008606.64947441
( L6 J& v( k' t& r0.735  秒

% ]) Q2 D- m! h/ t8 j8 u可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

2 ^. E! o1 |6 b2 E－－－－－－－－
1 D: |6 v# y$ P8 F9 ?1 w% V
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。