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

forcal

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

1. //用C++代码描述为：
   ) f, a: R9 m- g. j
2. s=0.0;
   , m0 l+ u# D# L* @, P4 Q4 Q8 Q+ {
3. for(x=0.0;x<=1.0;x=x+0.0011) # B& q6 a( I6 U% q9 M2 ~, `* u
   
4. {
   + g' F7 Z# _, ]6 K7 ^
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   / |, a# Q8 l' N; j; f4 k
6.    {
   ) O7 B( K4 `' t
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));/ \' a\" ^* N( ^1 ]5 H: v% c8 V
   
8.    }
   $ B; z( e7 g6 M: F* x2 X) w+ E7 I
9. }

复制代码

Matlab代码：

1. tic' w3 T\" \7 ]4 G% ?7 e# V$ M' U
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   * O  O) h2 P: B) K1 y4 s* l0 h, S, j$ u
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ; u) W0 f) ?  E! e9 D
4. toc
   4 T: H0 ~, E) [% P; b) B( q) @
5. : g& K8 Y! I' e  _9 r# S' B* z, m
   
6. s =
   1 @# ]) L# |8 ^& D' J1 Q4 @
7. 4 S# Q- ]  d- X, g
   
8.   1.0086e+006
   6 w# r3 M/ @% r% y4 j' p' ?4 r& E
9. / P6 Q- v: J\" j: C- g
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];& ~! o, R: D$ D; [, v2 l  _8 @
   
4. mvar:& H2 ]\\" n\\" C) U! d
   
6. t=clock(),% V! R0 S9 r' [( I9 F8 M2 i5 }
   
8. oo{/ A2 b1 \8 h3 e- w0 R
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],6 c# E: V2 `) \6 h- r$ Z6 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]- k# u7 S\\" ]: c& v+ D/ c& W\\" b/ R
    
14. };
    
15. # @/ B2 q  o# I6 i+ n7 Z$ q8 K
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">& ~! o, R: D$ D; [, v2 l  _8 @</font><br /> mvar:<font class=\\"jammer\\">& H2 ]\\" n\\" C) U! d</font><br /> t=clock(),<font class=\\"jammer\\">% V! R0 S9 r' [( I9 F8 M2 i5 }</font><br /> oo{<font class=\\"jammer\\">/ A2 b1 \8 h3 e- w0 R</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">6 c# E: V2 `) \6 h- r$ Z6 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\\">- k# u7 S\\" ]: c& v+ D/ c& W\\" b/ R</font><br /> };<br /> <span style=\\"display:none\\"># @/ B2 q  o# I6 i+ n7 Z$ q8 K</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641

9 p' [, Q; }- {  h3 b9 KForcal比Matlab稍慢些。

) E- e! m  H% X% Z: U3 h  ~4 O－－－－－－－－－－
. q" _6 C5 u7 p! D, D# ]
再看循环效率。
' k; j: x5 z, i: u" N' ^  s+ {
3 V6 x6 U" n! ]: Z2 WMatlab代码：

1. tic$ e\" ]3 @! S% T: U3 m
   
2. s=0;
   4 L, A' o, M' y3 i# I3 @
3. x = 0;* n2 y( A+ F) w' z\" k6 f2 K, h6 L) {5 M3 C
   
4. y = 1;4 l+ [% k; Z0 O; ^
   
5. while x<1    3 s1 M& J6 G- K2 X/ B
   
6.     while y<2;        , ^- u3 f& x- I
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 I& r; m; b\" Z$ M$ e
8.         y = y+0.0009;        
   5 ?  X, y\" ~! O; e
9.     end
   . M/ A+ I' U$ H& H' B0 M
10.     x = x+0.0011;
    ! C1 b% h: x4 r* N0 n
11.     y = 1;
    0 u& H$ I4 C, S- w: }5 k/ C4 l
12. end
    $ m5 E* Q8 [0 s& ?: _4 R
13. s
    1 s8 P; N/ B3 ^2 }4 T7 W+ m+ u  [
14. toc
    ! v  ]) p- y* M) s2 Q3 v
15. 
    ' p/ g) L# [- H1 `' O1 ?9 k
16. s =
    ! b3 r5 F8 d7 b. ?; f. y' v\" G
17. 
    8 D* c1 Z7 U: Q0 m# @
18.   1.0086e+006* @% W; K+ C\" v5 ~) V3 j7 @
    
19. 
    ' k5 g\" E3 g& G
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:2 p4 ^' F. Y  f' P/ o
   
2. t=sys::clock();
   ) z4 V8 O- G  C1 h& K
3. s=0,x=0,
   . R. e7 t: B! u$ T* E, ]\" U
4. while{x<=1, //while循环算法；
   6 G+ U) c; E1 U) O) ~6 b
5.     y=1,
   0 k1 B4 b! n8 s! b0 ~% ~$ x
6.     while{y<=2, 2 U0 \: R; [7 }
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ) ?+ c5 w. l8 D: G, v
8.         y=y+0.0009 8 f* |0 T3 c( J# u. N* A0 ~
   
9.     },
   + J4 S5 t0 ]( T\" z; M3 ?9 I6 Q7 g
10.     x=x+0.0011
    # s* F% B\" x) Q# h8 l' E% T
11. }, 2 I* m6 O. \1 S, z( S
    
12. s;1 j6 J9 U% W\" d1 R
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
" f1 a( |2 D8 x2 J0.734   //时间，秒

& E- o- [" s# Y/ b7 v我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
+ L6 _1 @1 B& `2 j
－－－－－－－
' D0 R& G9 @3 ?, A- E
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:  z8 e, P) W9 h& h% o
   
2. t=sys::clock();( C8 J$ M* m! R* t- q2 C7 M
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); , x, J0 }& _3 |& I) |! w
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   0 n. A\" V( U4 Y7 s' [/ F\" H
5. [sys::clock()-t]/1000;

复制代码

结果：
3 R. I" ~) k8 h1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
9 G0 l( F/ n/ [# v
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   + m9 j5 u1 ]% C: _3 X4 ]
2. tic
   * r: R# A  I4 p: S& v: r
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   & l2 z$ |( Y/ k, f
4. sum(sum(arrayfun(f,x,y)))( A4 N: w9 i5 q9 r6 D0 D
   
5. toc3 I6 v5 p; i+ `3 e9 Z5 b
   
6. 
   ! P& P. O' o6 X* V
7. ans =+ I$ ]0 F( @% U6 v% Q
   
8. ) t3 Y8 s8 c7 J: o4 A4 Y
   
9.   1.0086e+006
   1 L( Y  A2 H( j& E9 Z
10. 
    ; [* H, ~3 y! E, A$ V) G
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. % `0 |  m& Z- B4 R0 N% Z5 D: W5 i4 K
4. mvar:
   
5. ' K; e5 D# k: w4 `/ Z\\" r# b
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 3 W3 p6 t2 r9 L\\" u  r6 Z4 I
   
8. t=clock(),
   
9. 0 W- j# J: T3 q\\" ~6 S1 |  K6 s
10. oo{\\" B$ z- h9 d\\" ~/ @+ l
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. / r( |! _8 Z- ^9 s
14.     arrayfun[HFor("f"),x,y].Sum[]3 o+ ^4 o7 o6 \$ r8 E4 a
    
16. };
    
17. 1 W+ B+ N4 L\\" n- q  E2 ]' m# \3 c
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">% `0 |  m& Z- B4 R0 N% Z5 D: W5 i4 K</span> mvar:<br /> <span style=\\"display:none\\">' K; e5 D# k: w4 `/ Z\\" r# b</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">3 W3 p6 t2 r9 L\\" u  r6 Z4 I</font><br /> t=clock(),<br /> <span style=\\"display:none\\">0 W- j# J: T3 q\\" ~6 S1 |  K6 s</span> oo{<font class=\\"jammer\\">\\" B$ z- h9 d\\" ~/ @+ l</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">/ r( |! _8 Z- ^9 s</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">3 o+ ^4 o7 o6 \$ r8 E4 a</font><br /> };<br /> <span style=\\"display:none\\">1 W+ B+ N4 L\\" n- q  E2 ]' m# \3 c</span> [clock()-t]/1000;

结果：
0 {. l; C. ^5 e1008606.64947441
% t( C, P. E8 @8 }' \- u& h6 m0.735  秒
8 J" q! K  Y: f& X
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

& ~# A/ A  u6 F' V, k从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。