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

forcal

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

1. //用C++代码描述为：( E0 J/ N9 ?) I9 [$ F
   
2. s=0.0; . S\" G4 m  f( i+ q1 L
   
3. for(x=0.0;x<=1.0;x=x+0.0011) * U! G% r1 L5 {5 ]& w
   
4. {
   1 O- `& ~% W0 R* J& m
5.    for(y=1.0;y<=2.0;y=y+0.0009)+ ]/ S; ?* }% ^7 B( M8 M% c
   
6.    {* t1 v4 K: f* k. A$ F2 R
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));& B' u* v) k9 I* {
   
8.    }: U7 {# [* i! Y8 f  f2 A/ _\" S
   
9. }

复制代码

Matlab代码：

1. tic5 q. N7 F1 h9 t  ~* A! i
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);$ |4 P0 a5 j\" W; ?* {8 L0 [9 k
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   6 }, D; D\" h' D
4. toc
   % m6 M9 p0 ~* M/ H4 u1 K
5. 
   ' g9 @# O  d$ L) D
6. s =) l5 u+ Q( y, g8 [5 E
   
7. 
   + o% L9 K, x9 u8 W+ A- q+ H! {+ w
8.   1.0086e+006) g; C& a4 W$ z8 Z
   
9. 
   $ q, e\" t+ h# S/ v# @8 C9 W
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];) Q% X5 `+ P# z9 P
   
4. mvar:
   
5. 9 G% s7 v+ ]4 T3 ^
6. t=clock(),
   
7. 6 C7 l' \$ ^' n  k& u  G) ^5 B\\" ?8 k
8. oo{. v% I& I$ s( h0 a) W+ ^
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],8 c- S5 F9 b6 x! a: ]) Q
    
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. 3 T: b) z2 O, ~: d5 _7 h# K
14. };
    
15. ! k; J7 Z  ~4 X7 U
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">) Q% X5 `+ P# z9 P</font><br /> mvar:<br /> <span style=\\"display:none\\">9 G% s7 v+ ]4 T3 ^</span> t=clock(),<br /> <span style=\\"display:none\\">6 C7 l' \$ ^' n  k& u  G) ^5 B\\" ?8 k</span> oo{<font class=\\"jammer\\">. v% I& I$ s( h0 a) W+ ^</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">8 c- S5 F9 b6 x! a: ]) Q</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\\">3 T: b) z2 O, ~: d5 _7 h# K</span> };<br /> <span style=\\"display:none\\">! k; J7 Z  ~4 X7 U</span> [clock()-t]/1000;

结果：
& C2 j' s  z7 J) ?. H; U1008606.64947441
  g  P2 `+ g# z1 a1 J9 ?( R0.641
9 ?& @6 [, r! o3 ?4 M
Forcal比Matlab稍慢些。
' w6 p* Y1 H3 F
－－－－－－－－－－

再看循环效率。
, ^: o( e2 i' a" i1 ], S
* o, L9 e; h  p" n9 {( kMatlab代码：

1. tic
   ( I1 \7 U9 O5 a' ]+ }7 Y\" q) I
2. s=0;
   4 \/ Z  ?$ ]0 o; F# v9 L7 b\" l
3. x = 0;
   $ x1 x: @: W+ K* r4 I# N& m
4. y = 1;9 r! e8 Y) r) y. \5 W. |/ T
   
5. while x<1   
   2 a4 m3 P* w: m3 t: ~
6.     while y<2;        - d; e. j  C  G9 r8 J
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));/ @5 q) r' }- H9 R9 K
   
8.         y = y+0.0009;        / H8 F* b3 s. I) U' y, f
   
9.     end& \' f) k. K4 p% e) a
   
10.     x = x+0.0011;6 e& U- G' K7 U5 [\" e
    
11.     y = 1;' Z2 M- U& @7 _6 r. ]7 I
    
12. end! C1 u, `3 ]( s4 S7 M  @5 s0 }, A
    
13. s
    $ q: N( @8 f2 r- ?6 i\" B+ S8 Y( d
14. toc: Z; e8 O( `- k* G
    
15. 
    0 N4 m8 k) m2 Z' O
16. s =# q/ s6 a# c9 E; U5 h+ \
    
17. + U/ Z8 `4 V0 W; S
    
18.   1.0086e+0067 p3 _: N- G6 T0 ~; H\" m
    
19. 
    4 I; l( [4 ~9 g  A  G/ L
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:2 \, D0 q  t8 P& Q4 k/ a; D
   
2. t=sys::clock();9 D; J9 P' O3 g9 P8 C
   
3. s=0,x=0, 0 z. W1 T7 f5 i2 q6 G) w
   
4. while{x<=1, //while循环算法；
   \" P! s) V5 Q$ V3 r, k% J
5.     y=1,
   % b% I/ k3 C5 r; p; A4 j9 q3 @2 h5 ^
6.     while{y<=2,
     E2 M( j' a: F% f6 o& N, p3 H- o9 C
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), : |% a( p5 U* A
   
8.         y=y+0.0009 # g( `2 t8 C  c! g( E$ \
   
9.     }, 3 s' Q. Y# y1 M3 O2 T# y/ ]5 W5 ]
   
10.     x=x+0.0011 - L6 P# A) _9 L  q3 @% S3 i% n
    
11. },
    3 f5 {8 P/ Z: F3 i
12. s;
    % l\" Z1 c, R9 a( I* F5 E0 a$ I+ h5 I
13. [sys::clock()-t]/1000;

复制代码

结果：
+ x# n. M, d, T! n$ _. K; b' C6 @1008606.64947441
( j, M4 l6 O4 @7 |- P0.734   //时间，秒

& h) [' N, z% g/ ?" D* D9 }( M我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
: f& g2 \/ J" y( E2 [) Z. z. o
－－－－－－－
. f( \) }( ?6 A: Z: X" L& ~
& p7 Z2 P6 q' z# n  nForcal中还有一个函数sum专门进行这种计算：

1. mvar:; j8 [, m; e! |) z& H* E
   
2. t=sys::clock();
   : l9 X: n4 v1 n5 m
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ) |, g& _+ L4 m6 p
4. sum["f",0,1,0.0011 : 1,2,0.0009];* J' j- n. H. r. I\" b: T% |
   
5. [sys::clock()-t]/1000;

复制代码

结果：
. R' N2 f- \3 d1 O$ m1008606.64947441
0.719   //时间，秒

forcal

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

, t  s; G# {8 J; l3 [1 Zmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 C! a! i9 Q0 r5 d3 Z- d
2. tic: s' n( r$ y% ?  y6 ]
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);\" x. q+ l* A7 t. E\" |; N
   
4. sum(sum(arrayfun(f,x,y)))
   8 z# \% b2 b# }% t! F* }' M
5. toc
   . W8 B0 H- a4 h9 z* \% ^; b
6. ! x$ q- F1 ?- y8 z- _( @( {
   
7. ans =2 x' e, `' j* I) R- ?1 I
   
8. 
   ! b4 C) G3 [1 x; b
9.   1.0086e+0065 u1 _4 A6 o  k8 E' _, Q
   
10. 2 K2 i& m) ~\" b+ k
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 7 C, c+ F! X' r! H( \$ H
4. mvar:  L3 l; W  ?8 S# s- o- b
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. / M9 ]$ z+ j# `# S2 [1 k
8. t=clock(),
   
9. ; X3 ~! C\\" b! v9 e# A. O+ X
10. oo{+ R& n' t& N- \4 T, h
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. / J. H3 }, v6 t7 W( O
14.     arrayfun[HFor("f"),x,y].Sum[]6 `$ \3 N0 p% s) j+ d' z+ g6 m
    
16. };. Z# K- i; X/ Z: M
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">7 C, c+ F! X' r! H( \$ H</span> mvar:<font class=\\"jammer\\">  L3 l; W  ?8 S# s- o- b</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\\">/ M9 ]$ z+ j# `# S2 [1 k</span> t=clock(),<br /> <span style=\\"display:none\\">; X3 ~! C\\" b! v9 e# A. O+ X</span> oo{<font class=\\"jammer\\">+ R& n' t& N- \4 T, h</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">/ J. H3 }, v6 t7 W( O</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">6 `$ \3 N0 p% s) j+ d' z+ g6 m</font><br /> };<font class=\\"jammer\\">. Z# K- i; X/ Z: M</font><br /> [clock()-t]/1000;

结果：
! h: Y/ l4 p0 i1008606.64947441
0.735  秒

* ^/ J. M9 a# B9 H8 n# |% A# v可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
) ?; \7 Z8 ^: Y
+ e5 q1 `0 c( b* x2 E, I－－－－－－－－
+ D2 i* ?7 m3 [2 I
$ M6 Q$ S* Y0 D  @9 @从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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