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

forcal

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

1. //用C++代码描述为：/ E8 x0 Y  Q; V( F! D/ |$ S& J
   
2. s=0.0;   _$ |) D\" P( I0 r
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   1 L# F8 p, @7 u$ L3 r
4. {
   7 k: e5 G) T: L4 P
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   6 d; F: m) p0 O& `! ]
6.    {, L( x3 Z; l! \
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & y) @/ H2 g4 v! g5 N' F+ m& p
8.    }; W/ T  b\" K2 o8 n, g5 W
   
9. }

复制代码

Matlab代码：

1. tic
   5 Z, \- Y\" ~1 E) }7 U7 l9 `6 S
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);1 V. ~  r3 h9 u+ ^4 R
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))1 E; @; q: m3 Q. t1 q5 {
   
4. toc* Y/ ^& h. c3 [
   
5. ) h+ A! O; i. b0 h' y& _
   
6. s =
   * ]* t7 @0 u- q' e- \
7. 
   3 z2 o% _; O! B- k: d9 S9 `6 G
8.   1.0086e+006
   2 z! G6 H3 R4 T( g+ l4 f5 m\" n: _
9. 
   0 F. ^& C  l/ a# X) {$ Z& I
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];, {% J% ~0 g! N' q& g! n3 n: a
   
4. mvar:, o4 i! O# v0 F! f
   
6. t=clock(),8 a2 ~  `9 Q) A/ e4 T3 D2 e
   
8. oo{
   
9. * S6 _: v- o8 I# h
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],+ Q\\" q5 ?$ r8 ], b; x
    
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]; D% k: ]0 y) p# Z! C
    
14. };
    
15. # d. R' `; J; E9 l( R& `8 ?
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">, {% J% ~0 g! N' q& g! n3 n: a</font><br /> mvar:<font class=\\"jammer\\">, o4 i! O# v0 F! f</font><br /> t=clock(),<font class=\\"jammer\\">8 a2 ~  `9 Q) A/ e4 T3 D2 e</font><br /> oo{<br /> <span style=\\"display:none\\">* S6 _: v- o8 I# h</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">+ Q\\" q5 ?$ r8 ], b; x</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\\">; D% k: ]0 y) p# Z! C</font><br /> };<br /> <span style=\\"display:none\\"># d. R' `; J; E9 l( R& `8 ?</span> [clock()-t]/1000;

结果：
: O( C4 X2 O7 L4 z- N1008606.64947441
0.641

Forcal比Matlab稍慢些。
. {0 P  A2 h7 G3 D
- @% I( S5 B! @1 q& v  r8 D  I－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic
   \" Q  T9 ?/ Y+ m7 y, i0 M+ I) I' e
2. s=0;
   7 ~4 n# j) g- j  Y2 V# \3 t6 o% ]
3. x = 0;0 R. T- }3 x1 J
   
4. y = 1;( j' O5 p# b6 f) m\" X3 n8 R2 s& O+ I- l
   
5. while x<1    % U8 P0 N, x& S4 r\" }  b- a  {4 ?' m
   
6.     while y<2;        
   $ ^7 \, x( Y: G' S) b. m
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));6 |& j% c' i' ~( F; K+ ]7 C: o
   
8.         y = y+0.0009;        
     G- ?  Z! C: ]  e7 h8 F
9.     end, O- R\" h8 j7 `
   
10.     x = x+0.0011;
    ) H: W% M8 N1 o+ v$ o/ g! T
11.     y = 1;
    1 T( p0 s+ h9 U$ F  M3 C0 l0 R
12. end8 _$ D/ f, B7 q) _\" }
    
13. s+ s$ B% L, J1 |6 a% a3 O
    
14. toc
    ! ~6 U4 k! E$ B4 X  n
15. ; W$ I% H7 B& P+ B
    
16. s =
    . ]2 C/ s% B8 d) L9 M5 z
17. 
    - a+ @# ?1 O5 |+ h- ~
18.   1.0086e+006; K) n% k0 F& w* @' I$ h$ z0 l: e2 [
    
19. 
    1 H6 V6 s& ]8 T8 w! g! A- L
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:' a1 i  w5 @+ ^
   
2. t=sys::clock();2 N$ G$ C, q7 O! j3 O0 j  v
   
3. s=0,x=0,
   + K/ n$ w* r% [- z; H! U
4. while{x<=1, //while循环算法；
   2 [. X3 N3 k: o! X4 k
5.     y=1, : c( X$ s- ~\" U/ D( f2 O% |
   
6.     while{y<=2, 4 Z# |0 J9 V1 ^# L\" d: O
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 4 O1 b- ^\" a- }' n
   
8.         y=y+0.0009 ( i8 \+ @+ Y0 O( [* {& N3 Y
   
9.     },
   ! M! U& d' F% l. Y$ d$ t0 [
10.     x=x+0.0011 ) O* _, C7 d$ q; P4 q
    
11. }, : U1 I2 ^& t! [/ E& l  {2 L
    
12. s;0 i# L7 @' x9 U; v# [\" q+ D3 S
    
13. [sys::clock()-t]/1000;

复制代码

结果：
. s& Y2 F- g& ?/ h# |- q1008606.64947441
0.734   //时间，秒
5 `1 L; K, \, G; _7 M2 ]/ C
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

3 \8 T5 y) @! Y3 q－－－－－－－

Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   * b! y+ K+ W\" B+ n6 v) r4 t4 W
2. t=sys::clock();/ o' J\" x; t: N0 _8 ^' B
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( c; T# _# }$ ~/ O: X
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ) P; M9 X2 a8 L5 b( F
5. [sys::clock()-t]/1000;

复制代码

结果：
/ M5 [* `, I; M1 j1008606.64947441
7 u9 ?8 B6 ]8 z: w, W$ z8 H0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));; q% Y, h3 P. `. N1 I
   
2. tic
   % z9 ]7 L7 L% e1 y1 D* C* v
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);+ F/ E4 c0 W3 l$ r4 }
   
4. sum(sum(arrayfun(f,x,y)))
   ( I. Z: u\" Y3 A2 {1 k\" c+ E
5. toc- Z\" e  H8 H$ V3 ?: y% b; Q
   
6. . s# v! D) a1 V. m3 w, o. H
   
7. ans =
   \" S8 |/ m+ _* W5 @# {
8. 
   6 E  k: d$ n: C9 f. w! L
9.   1.0086e+0068 p% K\" g, w4 n7 J7 }& L
   
10. ( x( t  v5 n9 N2 w% @: P
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. \\" q. R+ N% t3 j8 ^+ L
4. mvar:
   
5. % S/ s! K# m. j\\" h8 G
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); - ~+ j\\" F0 y7 ?; f. k9 m$ d
   
8. t=clock(),
   
9. * m8 L1 @- E9 p: q- [- q& {  u  T
10. oo{  x( n\\" I7 r2 ?& S# b9 j
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],; L# i4 x: Q! L7 b9 B/ V
    
14.     arrayfun[HFor("f"),x,y].Sum[]% Z0 q) r! D( B7 Q
    
16. };
    
17. 6 p. d6 b  [  @; }9 d/ k
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">\\" q. R+ N% t3 j8 ^+ L</span> mvar:<br /> <span style=\\"display:none\\">% S/ s! K# m. j\\" h8 G</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">- ~+ j\\" F0 y7 ?; f. k9 m$ d</font><br /> t=clock(),<br /> <span style=\\"display:none\\">* m8 L1 @- E9 p: q- [- q& {  u  T</span> oo{<font class=\\"jammer\\">  x( n\\" I7 r2 ?& S# b9 j</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">; L# i4 x: Q! L7 b9 B/ V</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">% Z0 q) r! D( B7 Q</font><br /> };<br /> <span style=\\"display:none\\">6 p. d6 b  [  @; }9 d/ k</span> [clock()-t]/1000;

结果：
, g% n2 C9 N% u9 z# }, V1008606.64947441
" P: F' j! g# b5 O. ^6 H( ~0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
! X2 W) K+ y. q% O5 t, f
－－－－－－－－

从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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