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

forcal

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

1. //用C++代码描述为：
   ; l; W( U  l& ]; A
2. s=0.0; , [5 G* Z; u\" Z& ^\" h( [- ]2 G
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   * X2 R2 a# Q8 @5 h# l& a! Y# ~
4. {4 m, K# m- \8 a0 G
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   # ?( ]5 Q1 t; c/ u$ M$ \
6.    {
   2 t6 T; i  i$ A  V
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 w& s\" j! b+ E5 \2 P
8.    }) [8 N( u1 J: }3 ~7 L
   
9. }

复制代码

Matlab代码：

1. tic- ]2 c7 |+ F  H& \
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);7 L2 f! `5 @6 f! E8 v# U5 a' E; {' y1 E
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))5 D7 f1 d! w  w, H' Y
   
4. toc
   \" a; D( {) K( r9 U2 J
5. 
   / h1 X7 i+ ?# S
6. s =% O. H5 C9 g# ~) l1 ^/ h3 Y* G* m
   
7. 
   0 P! q/ C9 L4 o4 P\" U
8.   1.0086e+006! u  b. Q6 ]* R\" x
   
9. 9 g7 [+ D8 A2 L\" F
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3.   a\\" t/ n( z# l
4. mvar:% V' C9 I% K* |. h( q3 U: D6 h
   
6. t=clock(),0 `2 A  D8 s* l! m* G# _\\" T# z
   
8. oo{
   
9. 2 l! t: |( [$ m! K0 P* G  ^
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 A' T1 u' k\\" f5 q3 y+ I
    
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. ) T; a- U1 b' Z) p' F6 P
14. };9 G' J; r* b* c
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">  a\\" t/ n( z# l</span> mvar:<font class=\\"jammer\\">% V' C9 I% K* |. h( q3 U: D6 h</font><br /> t=clock(),<font class=\\"jammer\\">0 `2 A  D8 s* l! m* G# _\\" T# z</font><br /> oo{<br /> <span style=\\"display:none\\">2 l! t: |( [$ m! K0 P* G  ^</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">1 A' T1 u' k\\" f5 q3 y+ I</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\\">) T; a- U1 b' Z) p' F6 P</span> };<font class=\\"jammer\\">9 G' J; r* b* c</font><br /> [clock()-t]/1000;

结果：
# _6 T+ s- `, O' F5 G$ B' v* a1008606.64947441
0.641
) }* k5 [5 q+ i
Forcal比Matlab稍慢些。

－－－－－－－－－－

) c8 C2 J; S! S. K5 O再看循环效率。
1 M6 G1 V  |* U2 t9 ^
Matlab代码：

1. tic
   ! Y0 D8 P$ ~% {) v% {; _: y3 ~
2. s=0;9 Y2 f& g2 M* C9 ^
   
3. x = 0;
     O, W2 m) i9 K  l3 P9 x
4. y = 1;! s' e9 ~; v5 w( }& b4 A7 Q- a. C
   
5. while x<1    5 s5 d& \! h9 o\" I' z0 ?1 _6 c' a
   
6.     while y<2;        \" q8 `$ Y4 M: R# O/ A; T( h  ~
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));2 r% J& I3 Q3 c+ O* a$ E/ ~
   
8.         y = y+0.0009;        
     \, p0 _5 g3 h/ w* s
9.     end3 i6 `! X% v, T+ O  I
   
10.     x = x+0.0011;  \+ u5 G- U& \' c' I2 e% Q0 p' f
    
11.     y = 1;
    $ P# s4 ?# G% I9 {9 F( M
12. end/ P- F' c+ `5 x0 j( r* p
    
13. s
    4 L( w( A* B1 Q6 _9 V
14. toc
    ; ?  ]8 {, o, C
15. 
    3 M* A0 e/ T/ J6 s$ {7 d( q9 ~
16. s =' s$ d. c  t' E/ ]& E/ z
    
17. 2 l8 f' O+ e7 T7 h: k$ D
    
18.   1.0086e+006- l5 m9 Y6 E6 ]8 C' t+ t
    
19. 
    . S5 i6 D/ I5 R0 a' G$ v# T* Y
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   4 J7 z  d- B1 @, y; Z2 Z2 P- U
2. t=sys::clock();
   1 O1 B$ M( S4 A7 _  e: C
3. s=0,x=0, : U- u; W7 M+ d! G7 P
   
4. while{x<=1, //while循环算法； . x\" F- A% p% x/ v
   
5.     y=1,
   0 m# R0 o) Y1 p
6.     while{y<=2,
   ( B/ h0 S8 E$ J6 H& d/ B
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ; b( R; ?: D) q4 b& b7 j& E* q0 T
   
8.         y=y+0.0009 $ b! \$ |7 J- p. l
   
9.     },
   + s% |( N/ Y2 N
10.     x=x+0.0011 5 N9 q* u: R* P: g
    
11. },
    & Q. I- @! u3 x5 |$ G0 y
12. s;
    7 N5 K7 s1 @* {* r
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
7 o; G5 n* J! {5 Z+ O' h
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

& G: k7 G( q/ G7 Q2 Y3 r# m－－－－－－－

  l4 F( T* h! j9 q& e, _  n4 C7 LForcal中还有一个函数sum专门进行这种计算：

1. mvar:- q5 x3 T4 J8 X0 g
   
2. t=sys::clock();& W( o; z. ]! H4 m
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 U8 R5 _4 G: J& b\" `4 X% P+ b
4. sum["f",0,1,0.0011 : 1,2,0.0009];/ p$ v. e; l* [0 o4 |) b/ e4 I
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
( c) \1 k% v. O. c+ d# e& T+ Q( u+ ~
+ ^' O8 n* W) ^4 ^matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 h: ?' D9 ]) M  r: O7 `. J- p: z, J
2. tic
   ' G  U, C. y- g( @' {1 G+ c
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   \" T3 o\" R4 l9 i0 |
4. sum(sum(arrayfun(f,x,y)))
   . R; F3 @9 s1 V% @3 i; _  I; w& o
5. toc6 t  z6 \3 q+ _  R8 n; Q% W
   
6. 4 Q( `\" i6 M' `: h
   
7. ans =
   9 o' G: K( G/ y6 c4 J+ J\" r
8. 
   2 `0 g2 s. |! E4 Q* B5 u
9.   1.0086e+0060 O$ `6 L! y( I' ?) h! K
   
10. 
    $ V( z/ l. X1 n2 b
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];  F/ z& i8 n+ D
   
4. mvar:
   
5.   z- R; u2 ~. |4 @/ t
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. , T1 n# G; f5 K4 v
8. t=clock(),
   
9. / v! ~; [2 Q0 B7 V. m
10. oo{
    
11. 9 |4 g5 v& L2 s+ O( _
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],6 W8 m/ z3 H! o) \- {( V
    
14.     arrayfun[HFor("f"),x,y].Sum[]+ s7 i1 [# P* f9 _9 T
    
16. };( U: ~& t6 N3 O. P
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">  F/ z& i8 n+ D</font><br /> mvar:<br /> <span style=\\"display:none\\">  z- R; u2 ~. |4 @/ t</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">, T1 n# G; f5 K4 v</span> t=clock(),<br /> <span style=\\"display:none\\">/ v! ~; [2 Q0 B7 V. m</span> oo{<br /> <span style=\\"display:none\\">9 |4 g5 v& L2 s+ O( _</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">6 W8 m/ z3 H! o) \- {( V</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">+ s7 i1 [# P* f9 _9 T</font><br /> };<font class=\\"jammer\\">( U: ~& t6 N3 O. P</font><br /> [clock()-t]/1000;

结果：
/ a8 o3 k$ w. w, H  s1008606.64947441
0.735  秒

- T' f0 R9 J" o# H! I# B可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
& E/ k! R( p' Y/ I# J
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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