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

forcal

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

1. //用C++代码描述为：
   1 g  |  a8 u0 k! O; o
2. s=0.0; , J; N) o1 ?/ H\" y9 X$ j
   
3. for(x=0.0;x<=1.0;x=x+0.0011) 8 @- j5 T1 k  c: o: i# c1 V3 o8 B
   
4. {
   % q# u  B1 x. S4 O% ~  B, G. d
5.    for(y=1.0;y<=2.0;y=y+0.0009)5 @' Q* p3 P\" q+ f6 S! l8 e1 o# {+ I# V
   
6.    {4 m- i& u) o( e: a1 \
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 q9 Q! G4 U+ r3 ]
8.    }% K9 }7 n; C1 U( s) L9 U
   
9. }

复制代码

Matlab代码：

1. tic
   , @# O( p  }! L+ k\" X
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ! b$ Z) C: k! p! W2 R' W/ Z2 G
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   3 |1 s2 U1 e, U+ E; _
4. toc
   \" @- A( m5 r* o6 a  U
5. # z( q2 \4 m& A+ e! [) e
   
6. s =1 B. L7 [# _\" p/ e) _/ T- P' m
   
7. . n8 O- p/ d% M4 l. y4 i; Y
   
8.   1.0086e+006
   9 u; M. w/ s6 {1 t8 A/ ]7 K7 Z7 q
9. 
   1 T) s) ]) ^* k4 R1 r% ~
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; O2 E9 @$ ]+ {3 A/ u7 C
4. mvar:
   
5. ! H& e' u/ L2 b$ }
6. t=clock(),
   
7. , d' m2 Z# K3 D* a$ t6 V
8. oo{7 j. s. H1 t$ O6 W
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],7 C2 _/ W, J8 S8 j/ h
    
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]& N  Y3 a( {7 u5 Z6 p& r6 P
    
14. };3 P6 m) d$ I# K% k5 N! J5 L
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; O2 E9 @$ ]+ {3 A/ u7 C</span> mvar:<br /> <span style=\\"display:none\\">! H& e' u/ L2 b$ }</span> t=clock(),<br /> <span style=\\"display:none\\">, d' m2 Z# K3 D* a$ t6 V</span> oo{<font class=\\"jammer\\">7 j. s. H1 t$ O6 W</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">7 C2 _/ W, J8 S8 j/ h</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\\">& N  Y3 a( {7 u5 Z6 p& r6 P</font><br /> };<font class=\\"jammer\\">3 P6 m) d$ I# K% k5 N! J5 L</font><br /> [clock()-t]/1000;

结果：
; v. R& L8 f- H6 ?! _+ S6 O1008606.64947441
1 }- g3 w# o3 l* g) L+ e0.641
4 K, d* y" m5 d' H9 x
+ F2 m1 e6 x. H" B* v: s3 j0 ?+ PForcal比Matlab稍慢些。
* u/ {* v4 R* i; b$ n- a! g8 i1 L1 S
－－－－－－－－－－
) C# l8 p2 ?& v( z. y
; I, p. _  N5 ^" e/ D; l再看循环效率。
/ k- E# p1 o, {. A. F  A9 a/ s
Matlab代码：

1. tic
   + O- E0 C) d% G7 v8 W1 c
2. s=0;+ d4 A1 M9 ]: C
   
3. x = 0;* W  I! J: z. O8 w
   
4. y = 1;
   9 p$ k  V% x3 }9 s4 v* |9 q- u
5. while x<1   
   , J\" ?5 {5 A1 g) [0 S2 q
6.     while y<2;        * `. |+ S$ [) S+ i
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  F7 H: {( q( i- |6 @
   
8.         y = y+0.0009;        4 j' c) e8 D  ?- X% a
   
9.     end
   \" x9 k# I4 y8 f8 K& H
10.     x = x+0.0011;
    ) u# e6 g8 a9 E, g
11.     y = 1;/ J' f$ v& u) R+ q3 l' N+ M
    
12. end! ^9 \. ~: r3 }. Z% r6 L, N; q' @
    
13. s4 Y' b( k& v! A7 g* c& e3 H: L
    
14. toc
    * i\" z7 ~' u. R
15. ( _2 t' U6 {) @. O4 l) Q  Y7 ?
    
16. s =! n9 [+ {3 X+ E# u0 K
    
17. 
    , `5 H  J9 R\" V0 x& l( W% v* \
18.   1.0086e+0067 q% C/ _4 O+ i6 x6 s- _) t
    
19. : T7 Y- s4 j) W! T* R1 n
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
     E1 m4 S( e& n5 Z: f$ l% J
2. t=sys::clock();
   , R2 x1 z: s3 X2 z8 P* Y
3. s=0,x=0,
   ; f- @3 b/ g( F! W- ?# X9 C% X, k
4. while{x<=1, //while循环算法； 9 L) y' e  e6 ?& U9 b
   
5.     y=1,
   * {1 q  T/ k2 Z% K1 T. ^
6.     while{y<=2,
   4 r6 [7 k+ A# D- U
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 2 I& k4 t( M+ A! t/ u8 S& q
   
8.         y=y+0.0009 / R. g) h8 ]$ l+ `& U* O
   
9.     },
   . Y' Q. D: R' |' D7 h
10.     x=x+0.0011
    $ c4 L' i2 n) M1 L5 Z5 Z
11. }, 9 c% |( p( ^; c. ]  K
    
12. s;
    4 F- _( V% H# @- _
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
& n( ~5 T( d7 j; b, s6 O. g* X
, M( r! O8 {; R6 N' c2 w我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

% n4 q# P  b# M/ q% r7 e+ p－－－－－－－
/ w, y1 D+ g# E8 \0 Q4 G# ]0 V1 |
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   6 q' ]: b\" v% W: g; ^3 l2 [) m8 S
2. t=sys::clock();6 K+ }! h3 I5 {/ A# w9 l6 @
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 7 O. x6 t9 ]- h\" c7 [  A
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];3 E/ J- Z( K: o, F
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
) |$ k2 i+ Z! N6 O/ B0.719   //时间，秒

forcal

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

/ H, y5 U) a- }* Lmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ C8 U$ v& y9 x2 R& \9 Z. E
   
2. tic
   ) m% ^- l0 G5 j+ \
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);+ ]( a. w! e( Q7 h' x
   
4. sum(sum(arrayfun(f,x,y)))
   8 q! p1 p1 E4 Z: V* \
5. toc
   , W6 B$ z) c1 x
6. 
   ' t0 A, R* G9 ~- W  b
7. ans =
   8 c\" r: [) i  u: a+ z
8. 
   ( ^& `& |6 `; g1 P! e; t
9.   1.0086e+0065 j9 N+ J: [; [\" \3 l( T* r2 D
   
10. 
    3 k& Q9 ?- ?- }, x
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; _! Q$ `$ M; N& I2 e1 L4 C% w
4. mvar:0 ^$ K0 Y; o: P) H
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 3 Z% s2 ]( x0 K  H; G0 s6 f! N+ S1 Q
8. t=clock(),+ j/ M\\" _: ]0 e! ^
   
10. oo{1 z$ V3 |5 q; H# ]) O7 p  _1 b
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. ; h! F$ h1 m% N5 z, r9 u) h
14.     arrayfun[HFor("f"),x,y].Sum[]/ [\\" T* d  C* y. e+ I
    
16. };
    
17. ( r6 O/ u# q& ]; H
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; _! Q$ `$ M; N& I2 e1 L4 C% w</span> mvar:<font class=\\"jammer\\">0 ^$ K0 Y; o: P) H</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">3 Z% s2 ]( x0 K  H; G0 s6 f! N+ S1 Q</span> t=clock(),<font class=\\"jammer\\">+ j/ M\\" _: ]0 e! ^</font><br /> oo{<font class=\\"jammer\\">1 z$ V3 |5 q; H# ]) O7 p  _1 b</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">; h! F$ h1 m% N5 z, r9 u) h</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">/ [\\" T* d  C* y. e+ I</font><br /> };<br /> <span style=\\"display:none\\">( r6 O/ u# q& ]; H</span> [clock()-t]/1000;

结果：
1008606.64947441
: d  G8 [& S/ `2 `0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
- G8 a- k6 t+ K3 E2 k! l
－－－－－－－－

3 X% |. E! |: K! w8 q0 A& I% }) L8 n/ i从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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