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

forcal

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

1. //用C++代码描述为：9 \/ J0 b* }5 p: f- Y1 ]$ u
   
2. s=0.0; 7 o7 ]# e7 E6 o4 L1 w. k0 G! u, z
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   . `& H# f& d* P
4. {' Z2 P* T1 ]: I: g$ x9 y
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   - x+ Z\" L% p8 Q0 Q/ B
6.    {1 ]6 u+ d' x5 `5 C, u3 l) ^
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; H' b\" y! \5 }+ c1 t9 F2 E
8.    }
   : V0 C  y8 h( t9 \! U3 _( K
9. }

复制代码

Matlab代码：

1. tic
   - u( a4 ]; ^2 l\" B/ f7 p
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);& v. I2 X7 C. \. c( Z. W2 J
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   # q  A% O\" |. o/ e5 g* J
4. toc
   2 C) S2 E. ^' J2 y3 N1 K1 s2 P! ]
5. 
   2 B+ {( T5 O2 ~, g; {/ x' Z
6. s =
   6 G) m/ H( V+ i\" N: t; [2 d
7. 0 Q4 d1 \1 s( E8 r! i8 e
   
8.   1.0086e+006& N2 P$ g) D2 _\" B7 ~7 m: y& B
   
9. 8 S, Y2 N7 p- w7 f% ~+ `6 R5 c
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];. {3 c; v2 e1 L* P# D6 H
   
4. mvar:2 u+ ]! G6 q4 t$ |6 [! s
   
6. t=clock(),
   
7. 4 U! ~9 ]% t) d9 s/ r3 p
8. oo{) O) F  r0 r3 N, b; t! z
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],, d$ x4 n  b\\" L7 U+ w4 A* T/ E% S
    
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]/ D5 E1 n# R. B# t+ s5 W6 E, A
    
14. };
    
15. 5 Z/ ?8 L) ]2 p9 S! m
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">. {3 c; v2 e1 L* P# D6 H</font><br /> mvar:<font class=\\"jammer\\">2 u+ ]! G6 q4 t$ |6 [! s</font><br /> t=clock(),<br /> <span style=\\"display:none\\">4 U! ~9 ]% t) d9 s/ r3 p</span> oo{<font class=\\"jammer\\">) O) F  r0 r3 N, b; t! z</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">, d$ x4 n  b\\" L7 U+ w4 A* T/ E% S</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\\">/ D5 E1 n# R. B# t+ s5 W6 E, A</font><br /> };<br /> <span style=\\"display:none\\">5 Z/ ?8 L) ]2 p9 S! m</span> [clock()-t]/1000;

结果：
2 T! s' I4 \8 a; j# z1008606.64947441
, K; c" p) E+ b0.641

Forcal比Matlab稍慢些。
5 ?" b% m2 O4 C2 Z4 V  c
－－－－－－－－－－

再看循环效率。

1 s1 P$ ~" ]8 C. J& k: PMatlab代码：

1. tic6 d4 t' d* {8 T7 c% I
   
2. s=0;
   + `  d& c6 t$ I: v& d9 W' @9 h
3. x = 0;* P4 F; F7 R& B
   
4. y = 1;- H2 B, ~5 z! x
   
5. while x<1   
   7 O# q( I; l/ j
6.     while y<2;        
   - a1 ?% X& u% L( Y
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! l$ \- ?\" r6 h& K& N8 q' M* {
   
8.         y = y+0.0009;        8 ]+ j. V) N% l\" r4 ^2 _7 M$ ?& z
   
9.     end
   6 \; V4 z. C, k2 ?5 j3 ~\" y
10.     x = x+0.0011;
    - Y6 V7 Y2 y( W0 }# G
11.     y = 1;\" ~/ {1 Z, K2 H5 v% A
    
12. end
    4 _\" a6 u' ?$ N- y' Z3 S' z
13. s
    - C. h: r8 B+ T! y. _, n/ ]: }
14. toc4 U5 S3 q$ F% v) `4 n+ q  c! I
    
15. 
    6 v$ F* U# |. L! C. W7 B  x
16. s =7 X% ?4 D0 Q! M' t) T6 b, c, i' N
    
17. 5 a' x- m7 f* f4 g
    
18.   1.0086e+006
    1 n8 z% s/ A# F2 I1 v; j\" l- J3 L! M0 U
19. + O0 X& L4 [9 A6 a' Q7 @
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:) J/ C) U) A8 [( w2 \* L: ?) l6 R- z
   
2. t=sys::clock();
   # d% P* V4 r+ V  F! \\" Z% y\" k\" P& r
3. s=0,x=0,
   - H- y4 V7 Q7 c2 ]
4. while{x<=1, //while循环算法；
   % ]6 O* B: w! M% K5 z
5.     y=1,
   * u3 N  ]) H+ W( O/ s2 G
6.     while{y<=2,
   % h& T& i! I3 U+ x. J% ?\" ^6 K
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), % B  Y  `+ g% O6 H' U* K
   
8.         y=y+0.0009 3 [0 h- l3 P* @  w& L$ I
   
9.     }, ! r: f) Y: \& L, C$ t- n; O6 ~, b/ M
   
10.     x=x+0.0011
    7 R0 A* q# ^( X6 k3 M
11. }, + v+ b) ^; j, p
    
12. s;- a4 z; `. w- D\" a2 z
    
13. [sys::clock()-t]/1000;

复制代码

结果：
4 ?- G1 R! o5 b; s1008606.64947441
0.734   //时间，秒
9 @: P1 {+ K  E4 d( Z' C1 N
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

. h1 j6 b* \7 x－－－－－－－
# U$ o/ h1 Q1 k/ D3 F3 N
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   ! l. Z6 K7 O' }
2. t=sys::clock();4 Z% d5 B& h. F8 n1 P' f
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" \. M: R  x- l/ I2 `
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   3 J% z3 D! {' Q  M3 S( W# A
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

0 P- p! R  K1 t5 G7 |8 Lmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) ^( O0 p6 z6 g7 u
   
2. tic
   6 C9 k1 Y# t. H\" @. C4 k
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);& n& K' y  n  r$ B
   
4. sum(sum(arrayfun(f,x,y)))
   / [1 g, I5 s: ]
5. toc
   0 i8 o+ B7 Z9 c/ n  @; V
6. * |& u\" Q6 [# E3 r3 r6 b0 r
   
7. ans =) z3 j8 G6 ]: H* Q) }! v+ m! F\" ~. Z
   
8. 2 S: y7 O7 _; R- @
   
9.   1.0086e+006
   5 u6 h4 e: ?1 i- Z
10. 9 B, l$ a, e6 q9 Q+ |: @
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 1 ?' h+ K/ V( V+ Q. W& O
4. mvar:7 E- M6 v0 [' s0 ^5 d
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ( P' `- y. s8 E
   
8. t=clock(),5 i. p* j% R* ~* Q: W3 t
   
10. oo{2 Y% J7 ]4 V6 s
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 2 j9 ?. f2 G/ g& ^
14.     arrayfun[HFor("f"),x,y].Sum[]4 s* q: V2 ~( s9 `% `& f
    
16. };$ b$ C) U) j& ?: K' i- D
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">1 ?' h+ K/ V( V+ Q. W& O</span> mvar:<font class=\\"jammer\\">7 E- M6 v0 [' s0 ^5 d</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">( P' `- y. s8 E</font><br /> t=clock(),<font class=\\"jammer\\">5 i. p* j% R* ~* Q: W3 t</font><br /> oo{<font class=\\"jammer\\">2 Y% J7 ]4 V6 s</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">2 j9 ?. f2 G/ g& ^</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">4 s* q: V2 ~( s9 `% `& f</font><br /> };<font class=\\"jammer\\">$ b$ C) U) j& ?: K' i- D</font><br /> [clock()-t]/1000;

结果：
/ \. \4 F+ q. Z7 C1008606.64947441
7 B; i, ~6 v2 r  A# j3 a0.735  秒
& m7 ^, `/ b1 K# m3 f
% X3 @6 P  Y0 [1 h9 y可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
. D; |* ?7 Y$ T* N* U0 N
) C7 @! D6 |& t& ]) ?+ V－－－－－－－－
8 `6 }  P& |! c) z, `$ l. v
8 n" Z6 J- D# h# ~7 `3 G5 l从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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