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

forcal

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

1. //用C++代码描述为：( h% N: V1 m+ }- \\" a& H
   
2. s=0.0; % e, [0 R  T- F- _; u8 L4 \
   
3. for(x=0.0;x<=1.0;x=x+0.0011) ; i* o3 v4 B) Y' F\" ^+ a% V
   
4. {$ c; C& [3 F) E4 Y( z
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)3 i, S\" k4 D& s# [) L( @1 M0 ]
   
6.    {3 v. t5 e* t/ K1 k! l
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));4 e. Y1 D5 ]+ E( c
   
8.    }/ j% V9 E8 t9 ]7 f4 X7 B6 g
   
9. }

复制代码

Matlab代码：

1. tic6 [0 @. i1 l/ b/ E& Z% Y
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   & v8 x& Z( m  l8 {$ R- ?7 P8 m# }
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   4 c* s- O5 S6 A( j8 V! c
4. toc
   1 T7 A* R+ ]5 p9 x\" G
5. ! \+ k$ @( a7 x8 v
   
6. s =3 Y' c& J3 D1 \
   
7. + S( t' \7 i4 A
   
8.   1.0086e+006/ ^; @0 [* v; M4 g
   
9. 
   ( z7 t( |3 R. _6 m& o& p
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3.   F6 B2 V) v# H; T* [& Z2 J
4. mvar:6 O9 g! n% k& z\\" J* X5 G( _
   
6. t=clock(),# e  p3 u: m- v. L  ~1 n\\" ^
   
8. oo{
   
9. 6 i  ?! r* L6 C! E1 `\\" g
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: Q+ f# [* {( o# k- n% 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]5 }% B& _$ ]7 A& _
    
14. };
    
15. ) X6 L. B$ c' P- I* R
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">  F6 B2 V) v# H; T* [& Z2 J</span> mvar:<font class=\\"jammer\\">6 O9 g! n% k& z\\" J* X5 G( _</font><br /> t=clock(),<font class=\\"jammer\\"># e  p3 u: m- v. L  ~1 n\\" ^</font><br /> oo{<br /> <span style=\\"display:none\\">6 i  ?! r* L6 C! E1 `\\" g</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: Q+ f# [* {( o# k- n% 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\\">5 }% B& _$ ]7 A& _</font><br /> };<br /> <span style=\\"display:none\\">) X6 L. B$ c' P- I* R</span> [clock()-t]/1000;

结果：
1008606.64947441
9 v) Z/ g, x7 O0.641
1 g" Z5 W3 a& ^% y/ i
Forcal比Matlab稍慢些。

－－－－－－－－－－
8 Z+ G; Q. N! ?
' e2 N  V: h9 c再看循环效率。

Matlab代码：

1. tic  _$ C# K3 S& ?+ o) A8 i& p$ x
   
2. s=0;
   . C8 L3 [( H9 n! R
3. x = 0;
   0 E% U\" H/ q$ x+ ^& r
4. y = 1;
   ' n; y6 h: e; z7 N: |
5. while x<1    ; b- K( X+ ^$ ~# J+ k* |  k
   
6.     while y<2;        ) ]% \8 Y, T  n' h3 @
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ) v& x7 I# y9 i7 M+ r& d( e6 q! T0 ~
8.         y = y+0.0009;        & S2 n' e/ Z5 e+ ^% x
   
9.     end: n3 v1 A' e  T) m/ w
   
10.     x = x+0.0011;! x: {, [* f' @( p. {' H3 m\" l
    
11.     y = 1;. |8 |; V. q4 }) f6 L- i
    
12. end
    8 w2 _* Z, _0 K+ Z: R& `
13. s7 l4 k8 L) l: t
    
14. toc
    2 z$ O  L/ |/ e7 @, D/ I) K8 r# ^
15. ( S! Q4 b' @- H8 p9 b8 G
    
16. s =
    ) a8 d7 Z2 A0 ?- m! T/ p
17. % N! s9 h0 @; S: J4 v
    
18.   1.0086e+006
    ; _# \3 {8 ?& Z4 U
19. ) A1 Y7 d\" j- s) U8 W6 i
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:- ]  V1 X* D* U& j
   
2. t=sys::clock();
   ( g6 ]2 ^$ I2 i! K
3. s=0,x=0, * H- @  \. z* {4 }, s
   
4. while{x<=1, //while循环算法；
   ; E, J' {9 e% p9 W& O
5.     y=1, 6 Q( [0 J\" {/ X\" w6 j& ^  T) `2 a
   
6.     while{y<=2,
   6 j+ h' g: B. |( w
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   $ E+ Q. ~6 `. w$ O8 B
8.         y=y+0.0009 * D7 k, e7 B- r$ z6 r
   
9.     },
   1 b& Q% j5 Y# q- _/ C: P( N6 k
10.     x=x+0.0011 , l% V5 A% H\" q$ f. R5 @
    
11. },
    ! C- `; {# p2 K' N\" `! c& \( V
12. s;  N2 i! W9 d, ]/ l8 y4 |9 a8 f
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
* c4 A( o$ b. x1 O+ q2 _- k
) w: y9 s( V& j! x9 a, b& e我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
" X4 L6 o+ p. U6 [
－－－－－－－

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

1. mvar:
   8 X% w6 ?) I! i4 ?, h- I
2. t=sys::clock();
   ! r' f$ ]5 N. \/ T2 x
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); / u4 V$ q0 F4 }' t
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];, Q/ e  P% q! J( h% F2 C
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
+ y. a' ^5 a" o) u8 ~0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / U1 x! W; U* `
2. tic
   ) ^% s, c+ @: ~% L6 s0 f# g: [
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, v7 K  U, B: M) r: q  m
   
4. sum(sum(arrayfun(f,x,y)))
   7 R' C4 Z4 h: r9 W
5. toc! B  N0 F0 A. j3 }7 p
   
6.   t: T$ e* N7 @2 Q. w
   
7. ans =' J\" C& ^' C+ ?+ X6 R. x# ]! ^
   
8. 
   ( O\" h! |$ F% E. x
9.   1.0086e+006
   - n* s# E* W4 P$ J% ^; f1 \
10. 9 K- a9 K8 y6 v% w- s2 L# m
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];7 D+ o6 `; Z4 v; N  T
   
4. mvar:
   
5. ( d  @+ X\\" ]7 Z+ k) R& ~* ?
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. + G, h: R3 B- e4 p
8. t=clock(),
   
9. - u; @& [- e2 U, k6 k
10. oo{% _+ P) q% W) d$ D
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 8 @$ F: L* e% a
14.     arrayfun[HFor("f"),x,y].Sum[]$ D/ a1 w% B9 A6 ^
    
16. };. j5 n& ?6 W' d6 V
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">7 D+ o6 `; Z4 v; N  T</font><br /> mvar:<br /> <span style=\\"display:none\\">( d  @+ X\\" ]7 Z+ k) R& ~* ?</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\\">+ G, h: R3 B- e4 p</span> t=clock(),<br /> <span style=\\"display:none\\">- u; @& [- e2 U, k6 k</span> oo{<font class=\\"jammer\\">% _+ P) q% W) d$ D</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">8 @$ F: L* e% a</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">$ D/ a1 w% B9 A6 ^</font><br /> };<font class=\\"jammer\\">. j5 n& ?6 W' d6 V</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
& u2 x7 _3 t+ q1 Y0.735  秒

0 B. {  L: ~8 y1 Q5 t; \可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
7 J8 K, _4 g+ G% ~
－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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