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

forcal

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

1. //用C++代码描述为：7 t4 M, d' E  Q9 U$ W7 r
   
2. s=0.0;
   2 j' o0 {$ ^* I( {$ v  }0 |9 |
3. for(x=0.0;x<=1.0;x=x+0.0011)
   1 @5 U7 `, f- g* j) |. i
4. {
   ) h/ S% h( D3 B- A! K  X& N' Z
5.    for(y=1.0;y<=2.0;y=y+0.0009)# k  c& g* E& t) m& |
   
6.    {
   : I0 g; ]6 a4 j9 x
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) d+ M5 O0 \2 `% c' u
   
8.    }% Q0 V5 h) A\" O2 E
   
9. }

复制代码

Matlab代码：

1. tic. a: P  H' c  l; }3 n4 G
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ) ]0 D) [* o7 E& n7 S0 |
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))% Z& {$ A) Q9 T0 Q- ]( O+ j
   
4. toc8 }5 S9 O2 D  H
   
5. 
   ! F' w6 ^2 }6 @0 _9 ~0 d
6. s =
   \" A% e6 W/ P$ H; m7 Z! ~/ N
7. 8 j: j9 m& h8 M% c
   
8.   1.0086e+006
   ' P$ l$ f( u' U
9. 0 Y& C) C; T& A: }9 W, D- I2 k& _
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];0 q) J, d9 M. T
   
4. mvar:, h. b2 f\\" j0 c7 ^# X
   
6. t=clock(),
   
7. 0 D. [% L1 x5 d\\" @1 u
8. oo{$ ~3 E\\" X( Y  L# Y' \
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],' c: }: o) n/ J1 ^
    
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]\\" z, b+ T1 ]' |0 o( k* E, _9 F! d
    
14. };5 Q1 I  _, Z6 i/ S
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">0 q) J, d9 M. T</font><br /> mvar:<font class=\\"jammer\\">, h. b2 f\\" j0 c7 ^# X</font><br /> t=clock(),<br /> <span style=\\"display:none\\">0 D. [% L1 x5 d\\" @1 u</span> oo{<font class=\\"jammer\\">$ ~3 E\\" X( Y  L# Y' \</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">' c: }: o) n/ J1 ^</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]<font class=\\"jammer\\">\\" z, b+ T1 ]' |0 o( k* E, _9 F! d</font><br /> };<font class=\\"jammer\\">5 Q1 I  _, Z6 i/ S</font><br /> [clock()-t]/1000;

结果：
1 q+ a+ q* x# p2 t0 h0 t1008606.64947441
0.641

2 l5 q- t* h& h- sForcal比Matlab稍慢些。

% C: m% v2 Q$ S6 w－－－－－－－－－－

再看循环效率。
5 x: @+ K4 Q: m; H  W. g
5 Z! x0 d4 Y- L6 O8 J6 \Matlab代码：

1. tic
   5 w\" s$ M6 c. `6 ^+ }/ p
2. s=0;4 z+ Z& b3 a. {7 O) q
   
3. x = 0;1 M! `% q/ V: Q* G0 Y8 D$ v
   
4. y = 1;
   * P' _$ T0 F' ~% e) R
5. while x<1   
   \" v6 c$ n& a! B6 L6 t: E  k
6.     while y<2;        0 n) z: S6 _2 A2 X$ g
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));( Q  L% ~' [3 j- N- d
   
8.         y = y+0.0009;        # e9 ^0 `  Q9 v. ?. E* z9 S
   
9.     end
   : }* t6 b7 T4 f
10.     x = x+0.0011;
    % q* q\" H3 V8 O$ N
11.     y = 1;; x$ ?5 l# w! c% q) g* ?
    
12. end
    - D; C8 x& O7 U
13. s
    . T2 j0 f- a) J; L$ z$ M& g, J\" [
14. toc
    0 U- w0 {) P) T9 t6 {- i0 F7 R. _
15. \" X3 |; C2 o$ H5 a* h  U
    
16. s =
    $ l4 G, ]: \. U\" G. C* U  X
17. 
    / X) a5 T: q$ X) E! p% Q- E7 s
18.   1.0086e+006
    9 f. Q& ?! a\" s+ S
19. 
    - {5 r; [\" ]% \) z* p+ I$ ?( J
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   ! L* i' z/ X7 W) j8 b3 Q8 |# Z
2. t=sys::clock();+ w3 F; a% ]& h. @7 t
   
3. s=0,x=0, : o3 c+ {, R) b
   
4. while{x<=1, //while循环算法； & H6 W. u# k, R0 w3 ~$ h
   
5.     y=1,
   + v, E7 D+ I  j
6.     while{y<=2,
   5 I  R' B9 A; B9 [; n/ f% [
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), $ f; B% G8 w' W) D: L; H
   
8.         y=y+0.0009 ) H, W! p+ r- ?2 O8 K6 X\" Q- K; I
   
9.     },
   3 j0 w+ J+ z: H6 j; {
10.     x=x+0.0011
    - e/ v, Q& l, U1 D
11. },
    : j/ W% [) e) x) P
12. s;
    ! q8 J& P5 r\" @  ?5 \- n: W
13. [sys::clock()-t]/1000;

复制代码

结果：
; i8 a# `" c6 z$ H+ _( @+ `1008606.64947441
& ?& P  {' {/ x1 g/ }! Y0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

. F' }1 z6 C+ L9 r* \4 S9 s－－－－－－－

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

1. mvar:
   3 ?! Z+ b- U: D/ a
2. t=sys::clock();
   / \( i9 x$ m1 h0 S; B  U
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' Y4 N0 A+ q; O5 o) _\" V4 i
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   4 e% L, o9 A3 }) _' g
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
; u6 }" q& [( Q0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
2 i! W1 m4 k! h- k  s
) Z7 }* ~# w9 [5 p9 ]# @& V  `matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));\" X2 {; U8 y( p
   
2. tic+ v$ ^: f* {( c3 n5 \7 x
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   0 g6 ?9 l7 @& Z  J$ R' J
4. sum(sum(arrayfun(f,x,y)))\" t& B. b2 g: I$ U0 h& s
   
5. toc% S6 O8 q7 C. f( r6 m1 a1 ^
   
6. + x4 Y0 K8 `9 n: }6 L
   
7. ans =
   ) z& G' c* `& P9 H3 V! ^
8. ; \4 t; u+ F6 s. x' d
   
9.   1.0086e+006+ T, q$ O2 O5 q% ~3 y8 p, S- a
   
10. - h2 o8 j# L, h6 S% f5 W, K
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];$ X3 B- F) u; p
   
4. mvar:3 H) {& z4 B, p4 F/ ?5 @, _
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); % Z( z4 `9 i' Z1 T' v0 K; O5 q
   
8. t=clock(),& [$ S1 ^$ R+ Z- u% A
   
10. oo{3 P2 i) U: B/ A: D
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],, k* E7 n: R# H: o7 y) N
    
14.     arrayfun[HFor("f"),x,y].Sum[]6 q: P: R7 ?& K* I( _$ ~0 I
    
16. };4 T5 U1 W7 s8 l
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">$ X3 B- F) u; p</font><br /> mvar:<font class=\\"jammer\\">3 H) {& z4 B, p4 F/ ?5 @, _</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">% Z( z4 `9 i' Z1 T' v0 K; O5 q</font><br /> t=clock(),<font class=\\"jammer\\">& [$ S1 ^$ R+ Z- u% A</font><br /> oo{<font class=\\"jammer\\">3 P2 i) U: B/ A: D</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">, k* E7 n: R# H: o7 y) N</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">6 q: P: R7 ?& K* I( _$ ~0 I</font><br /> };<font class=\\"jammer\\">4 T5 U1 W7 s8 l</font><br /> [clock()-t]/1000;

结果：
! J" j2 i4 L- F: f4 i1008606.64947441
0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

" ]3 r& G$ x: }4 T9 u! |－－－－－－－－
+ v( ~0 B; X5 O
3 b4 u! _5 r4 i9 v, N- G; \从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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