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

forcal

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

1. //用C++代码描述为：
   ! f* k\" H/ Q; s/ l- K\" ?
2. s=0.0;
     ?# e$ z% w! p- B, C
3. for(x=0.0;x<=1.0;x=x+0.0011) 4 @  S0 A; Q, q6 @' c- @7 C  v
   
4. {5 K$ I5 g2 Z, y
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   + F/ L5 _, b, r; u5 z; q3 F) Q
6.    {
     c9 p8 H& b% N0 o
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));. \9 n- ^; @6 C8 o% m
   
8.    }, x+ W; n, o1 S' F& j
   
9. }

复制代码

Matlab代码：

1. tic
   3 e0 T7 Z( Z$ f( N  ]\" n% u
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   # O4 H3 P( ?3 t8 B
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))7 n3 @  K0 ?- I% d4 |  h- Z% F2 `
   
4. toc1 U; a0 r1 g7 i  U+ [( F
   
5. & o2 q. |5 I: A, Q  g# n3 f: q
   
6. s =
   - G. @6 o$ G# l6 E7 [, ^& n9 p
7. 
   , c% }; R$ a; r\" k+ `
8.   1.0086e+006
   5 W2 A6 q- s( _
9. ' q$ X! X# P% N' ]: x4 J6 ]3 L
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];4 }1 l1 E' O; n7 P4 _. b: R& Z
   
4. mvar:
   
5. / L+ a' h# `/ F! `: N% O2 r
6. t=clock(),
   
7. * M\\" Y4 \0 V' D) }2 a' i$ P0 z
8. oo{! c6 \/ T4 t$ M3 C! H1 P8 [
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. $ N1 e- J( f, S( D\\" Z0 b% O: q: U
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]' [  O( H8 H  Y: ~4 }
    
14. };& Z2 b; B& A; A; p( l
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">4 }1 l1 E' O; n7 P4 _. b: R& Z</font><br /> mvar:<br /> <span style=\\"display:none\\">/ L+ a' h# `/ F! `: N% O2 r</span> t=clock(),<br /> <span style=\\"display:none\\">* M\\" Y4 \0 V' D) }2 a' i$ P0 z</span> oo{<font class=\\"jammer\\">! c6 \/ T4 t$ M3 C! H1 P8 [</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">$ N1 e- J( f, S( D\\" Z0 b% O: q: U</span>     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\\">' [  O( H8 H  Y: ~4 }</font><br /> };<font class=\\"jammer\\">& Z2 b; B& A; A; p( l</font><br /> [clock()-t]/1000;

结果：
3 _: Z4 X2 ?; F# K9 v1008606.64947441
" m6 V! v0 L% v7 [: [0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－
" ?* Q, B! \" _1 Y% E
3 h4 m9 w+ D# R. b7 T再看循环效率。
( B' |: [+ G' v& Z! N
# M$ p! Y7 m2 V9 E1 _Matlab代码：

1. tic9 t% w$ `3 M\" _$ O% P& B  V/ x& J
   
2. s=0;8 Y: ^\" ^4 O. |* W% ~9 `
   
3. x = 0;
     E: m5 n- y1 `. q- v' Y
4. y = 1;
   \" V. W8 ^9 |9 X$ R& Q4 n/ q
5. while x<1   
   - B4 M$ p: F\" {& v
6.     while y<2;        
   6 K# }; h. T8 G& c7 B
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));$ y' D; a/ R$ U: D' g; L
   
8.         y = y+0.0009;        9 q: c- |6 Q( K' F3 W* E  d
   
9.     end/ T! A: a3 T! \. T) q# c
   
10.     x = x+0.0011;
    \" x8 |& ?* Q5 c9 w, ~
11.     y = 1;& ]9 b4 C8 e9 Y
    
12. end$ ], a( I/ h6 [% N% ]
    
13. s
    \" i8 S) Y# a! _1 C
14. toc; i9 H- m& w- l5 X. W1 r' M
    
15. 
    6 X9 C& r0 @' X# \# g5 g; L
16. s =
    . z: q) r9 }) Z: s: k5 l\" w
17. 
    $ D8 }# O4 F# l6 ]' d
18.   1.0086e+006& N  C/ ^% [\" v' ]/ w% O% Y
    
19. 
    \" {# g1 f0 D6 r\" P3 q
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   / L# M\" @. M3 Z+ y3 l
2. t=sys::clock();3 {2 ]1 y+ z7 B- x0 g+ [
   
3. s=0,x=0,
   3 R1 X\" w0 \2 y- }! }- I2 k8 _3 P
4. while{x<=1, //while循环算法； 8 ~% e3 D7 P3 i! K
   
5.     y=1, % A, g# C\" R% ~  y( e: i0 ]4 l4 v
   
6.     while{y<=2, 8 I9 ~  J\" ^& c2 t! j
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   \" D* r6 C\" F4 T
8.         y=y+0.0009 ' {( A/ |& y8 V+ N* Y( L6 E
   
9.     }, / i2 t6 Z7 w0 i0 G+ D
   
10.     x=x+0.0011
    , p$ `0 U\" z9 m
11. }, % O' v. t7 ^2 t+ V/ S9 Y
    
12. s;: Q7 v* ?9 W! d  D  Q3 R
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
; y; z. y7 |* N0.734   //时间，秒
. X0 L2 n; o5 o# V2 T7 V8 x& h
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
7 a' l, l( k5 Y) d: \/ z" K
6 P5 q1 _/ m( s& R- W% Q8 Z－－－－－－－

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

1. mvar:
   . X1 v) |5 e+ T
2. t=sys::clock();
   7 }% u% x: w+ e, `( A; q7 F
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 2 n9 N\" ]1 c, Z7 d5 I9 _
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ! D6 i1 P/ f& \# \3 G& V( g  b* r. ~2 a
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
) u# y! G+ x) {8 B- q3 E+ n
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
     k6 G. c# D# z$ e% h% e- H8 C- e4 K
2. tic
   * \& R5 @  u5 b- w6 ?
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 _* C  i( [: D1 B$ A
   
4. sum(sum(arrayfun(f,x,y)))
   & ]4 o6 K7 N, ~0 ?1 `! f! n
5. toc\" {: A6 {' a. E- q+ E
   
6. $ h$ i) n! p, Z* i
   
7. ans =
   1 j3 O& d\" z4 R4 C$ d0 X
8. $ k  V0 J1 D* d3 d; F1 T/ T# N
   
9.   1.0086e+006
   % j. }/ M: U4 ?8 U
10. 
    ( }' q/ K' j5 a4 ^
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];! w0 I2 p8 }% ]
   
4. mvar:0 H( V1 ^; u! `5 j( A$ N
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 8 v1 n9 L  u2 A& ?. m! w' H
   
8. t=clock(),5 @) K, J/ Q% ?' m
   
10. oo{
    
11. 0 j1 d4 h7 e8 t) a# f9 ~
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 3 h7 K) p' o! i\\" r. W( b9 k6 Z$ W
14.     arrayfun[HFor("f"),x,y].Sum[]- G8 h5 [+ D. G5 d# `1 ?# _: L
    
16. };( g! O; E% y6 F0 \1 B( u5 |6 s# d
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">! w0 I2 p8 }% ]</font><br /> mvar:<font class=\\"jammer\\">0 H( V1 ^; u! `5 j( A$ N</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">8 v1 n9 L  u2 A& ?. m! w' H</font><br /> t=clock(),<font class=\\"jammer\\">5 @) K, J/ Q% ?' m</font><br /> oo{<br /> <span style=\\"display:none\\">0 j1 d4 h7 e8 t) a# f9 ~</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">3 h7 K) p' o! i\\" r. W( b9 k6 Z$ W</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">- G8 h5 [+ D. G5 d# `1 ?# _: L</font><br /> };<font class=\\"jammer\\">( g! O; E% y6 F0 \1 B( u5 |6 s# d</font><br /> [clock()-t]/1000;

结果：
/ f  y; M" n( P. v) y1008606.64947441
0.735  秒
1 n: ]( I, h0 t" R* D# U4 Y, [
" N% w1 Q$ w) F+ e# U  u+ j# K可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
3 K, r  M9 f! @
－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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