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

forcal

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

1. //用C++代码描述为：
   8 A- w! I- Y' I! ^9 @! l1 \' {' T
2. s=0.0;
   * y) V) u; l, m/ \\" b$ k0 m5 P
3. for(x=0.0;x<=1.0;x=x+0.0011) 5 T' y* r' A7 U( S$ |! D. b
   
4. {# e, I' b9 g% k6 U
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)\" U) B: v- }& y\" O1 W9 L) g
   
6.    {1 Z7 g0 z$ f  s! o- l
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ }5 t- H' _/ z4 A+ y: F# m
   
8.    }
   . |\" n; U1 f6 l; k4 i
9. }

复制代码

Matlab代码：

1. tic. f# _\" @7 i7 U
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);5 Y6 y\" {' P, V3 D
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))3 k4 w3 ~6 R1 q* v% s6 |
   
4. toc
   & k0 i8 V\" W7 L* H4 B3 `
5. ; h! }7 }\" @. d. a/ a
   
6. s =0 I. y; e0 l# y  Q  F# g/ [\" \
   
7. 
   + a\" C. X, Y& O; L( ]% P
8.   1.0086e+006
   3 B\" |5 k' ?4 U; M* c
9. . z- ^\" z7 u' @9 A# x0 |\" d5 i
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];\\" j0 D: R2 M$ j
   
4. mvar:8 S7 M) ^. m9 k2 [
   
6. t=clock(),
   
7. 2 f/ K  W4 u$ s9 B+ D
8. oo{# {; v9 L( b7 @
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],& K! V+ |3 U! C$ O2 I  i# `4 W, t  {
    
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]7 s' T; _$ G& W6 z4 z
    
14. };
    
15. & n7 q! o7 M0 e& K+ B/ v* P
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">\\" j0 D: R2 M$ j</font><br /> mvar:<font class=\\"jammer\\">8 S7 M) ^. m9 k2 [</font><br /> t=clock(),<br /> <span style=\\"display:none\\">2 f/ K  W4 u$ s9 B+ D</span> oo{<font class=\\"jammer\\"># {; v9 L( b7 @</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">& K! V+ |3 U! C$ O2 I  i# `4 W, t  {</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\\">7 s' T; _$ G& W6 z4 z</font><br /> };<br /> <span style=\\"display:none\\">& n7 q! o7 M0 e& K+ B/ v* P</span> [clock()-t]/1000;

结果：
& t, a0 m/ ?$ d1 i( `1008606.64947441
0.641

Forcal比Matlab稍慢些。
/ J. A- c+ j0 U% q* G* U: s* }+ e
% s* d! s4 Q2 s! k1 V' Z－－－－－－－－－－
# Q$ }: T. ~3 F4 F
再看循环效率。
' L9 Q$ u) I) P6 N% H
Matlab代码：

1. tic  [- T' o4 X. s: L; A
   
2. s=0;( b& d, u: ~3 J7 J4 D
   
3. x = 0;
   ( }3 C: u2 j! E/ f. R9 J6 C, D
4. y = 1;
   % I1 J. U1 b' a' h( _8 j5 l6 \9 w
5. while x<1    ; h5 F4 C5 P, O  [) ]1 y/ ?( E
   
6.     while y<2;        2 |$ c4 d9 \0 y' b) d4 J
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ u5 k; ~5 O6 f0 W. x
   
8.         y = y+0.0009;        
   % D- y: ]1 j5 K1 B
9.     end
   4 [* h7 S+ ~. P( |+ n1 m
10.     x = x+0.0011;
    7 y+ S& O6 C1 K3 X
11.     y = 1;
    6 p0 r4 _- u0 C% ?' o6 `- p' {
12. end
    - w6 v8 d5 e9 P* d
13. s
    ( d. i! Z8 p' V3 v
14. toc5 r; F+ `7 H* L, ]$ B
    
15. 
    & k\" n2 J9 a- h9 ~3 l3 c: c
16. s =
    / O6 v5 _5 A6 |( ]: f* h
17. 
    / R* g3 i- }6 l& E# j
18.   1.0086e+006
    5 X+ T5 T, T2 y  Y\" I- C
19. 
    ( ^5 c0 _1 [, y( W: Y& a  v5 t; [
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   4 {\" B; |' R6 ]: m* X
2. t=sys::clock();; z. H9 z$ p$ m7 G. M
   
3. s=0,x=0, 4 C, e2 R0 G% d# M; ^2 @\" d; Z
   
4. while{x<=1, //while循环算法； & j8 K$ C% F- I: ~3 u2 ]2 c& v4 ]  o
   
5.     y=1,
   8 }8 f2 o+ t4 m+ H\" m
6.     while{y<=2, 2 Q- s( q# D% [  @& ?9 ?. T
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   . `/ @# x% p* l( b
8.         y=y+0.0009 % h4 B3 @; L0 `1 {
   
9.     }, ; s4 ^. C4 W1 B* C7 o
   
10.     x=x+0.0011 ! X1 a5 D5 C: m' `
    
11. },
    5 x5 ^! ]/ S) g$ R
12. s;
    ! b' a% g9 D' {! w( l! r
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
' |* G$ @' [/ M3 k; m0.734   //时间，秒
' _) T& w2 x9 `! @
  ~4 M* f, @; \5 z: L" p我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

. G6 |; ?* r5 Z4 V－－－－－－－
" m  u5 D3 e& E: q9 l& n$ J
$ x9 j2 U1 ?9 ^* _( p4 u- UForcal中还有一个函数sum专门进行这种计算：

1. mvar:, E: K8 q$ H2 N& [5 e5 x  K! m. s
   
2. t=sys::clock();9 g6 G' w  ]* m3 z0 {  M! ?
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); : C: q. j: N+ d' o
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];% a5 Y% `' z9 g
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
* f/ x- v) o+ Y* J0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
, {# \: p& S9 W$ T  _4 ~
2 c2 v9 q# g3 [* X/ K% Jmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));' P( p4 ?7 Y1 [9 A
   
2. tic
   : I. r# |. {7 {\" P0 Z/ M; Q
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);5 Y  `  ?6 q( I1 H* g3 c& o
   
4. sum(sum(arrayfun(f,x,y)))) `$ A5 I+ n: y
   
5. toc
   8 ~8 A$ \4 I1 q
6. 
   - O4 d& k7 p( [# {
7. ans =4 |/ V5 _- h$ g\" h- ]: D, x
   
8. 8 u2 }2 b2 E' e* D( g7 t- ~' B
   
9.   1.0086e+006
     x1 l0 b4 f\" v' w% {
10. 
    6 T$ ~' Q6 x3 K: \1 N. d2 f
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. # z0 x+ A\\" ]/ a+ i9 x1 z& k
4. mvar:) t2 X. `/ Y\\" [& [
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ' x, {! D7 L8 N- }
   
8. t=clock(),
   
9. & M7 e& q) O\\" x8 {1 \7 F
10. oo{! t$ |2 ?1 O! T2 ?; r  U$ x1 }! F2 C  w
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 9 \: J8 _8 N) K9 [  w7 z
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. ; \# f; ^4 j4 I; l
16. };* F8 D( {, D% {, {, [6 ?; Z  S
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\"># z0 x+ A\\" ]/ a+ i9 x1 z& k</span> mvar:<font class=\\"jammer\\">) t2 X. `/ Y\\" [& [</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">' x, {! D7 L8 N- }</font><br /> t=clock(),<br /> <span style=\\"display:none\\">& M7 e& q) O\\" x8 {1 \7 F</span> oo{<font class=\\"jammer\\">! t$ |2 ?1 O! T2 ?; r  U$ x1 }! F2 C  w</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">9 \: J8 _8 N) K9 [  w7 z</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">; \# f; ^4 j4 I; l</span> };<font class=\\"jammer\\">* F8 D( {, D% {, {, [6 ?; Z  S</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒

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

－－－－－－－－
9 a0 l" V( s' [% d: G; o5 P
8 ^4 u9 h8 f( X5 F- X/ }+ X7 Z! X8 p从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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