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

forcal

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

1. //用C++代码描述为：. e4 G* e* ?\" o. Q2 V( N9 t7 @
   
2. s=0.0;   R8 U) n, {  b$ l9 R
   
3. for(x=0.0;x<=1.0;x=x+0.0011) # x  L\" X% T\" V! s% s
   
4. {
   ) ~; C  d7 X6 U$ W4 @
5.    for(y=1.0;y<=2.0;y=y+0.0009)- v4 I# `$ Z/ u) Q0 G( e
   
6.    {
   9 v  I. i/ ]. o! i( ^: D9 |
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ s8 y& H7 c1 O$ q' d
   
8.    }
   ; q( W# r1 t( |
9. }

复制代码

Matlab代码：

1. tic
   % x& v7 {) S& y\" l/ t/ e6 N
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);2 D! i; @- z3 m# `% e# V/ U; {- b! L
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))). Z7 B1 l/ @7 Z; W& b$ x& T
   
4. toc
   + C\" u% K+ s0 w3 u9 u% u4 I& a) R
5. 
   , ~% I4 ^# `- R$ e
6. s =
   ! |7 B) Y- T4 q+ B\" x+ t' A' V4 F
7. 
   # |$ a0 i# N% Y. w& d9 S2 O( r
8.   1.0086e+006
   4 J* d0 {; E, F/ L, ^/ v& c, a
9. # q; Y  G$ M6 m& Z
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. / C7 `4 J4 K' `& O
4. mvar:
   
5. $ Y- w% i* S4 t1 Y
6. t=clock(),+ y# m* G0 g, |7 _( o  X/ B
   
8. oo{
   
9. # S+ i2 ]$ j! O) j  p
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. * i% }; ~3 ^/ h) n
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]6 }5 {2 v6 ^3 C. l
    
14. };
    
15. # W2 |0 e6 \) P+ r
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">/ C7 `4 J4 K' `& O</span> mvar:<br /> <span style=\\"display:none\\">$ Y- w% i* S4 t1 Y</span> t=clock(),<font class=\\"jammer\\">+ y# m* G0 g, |7 _( o  X/ B</font><br /> oo{<br /> <span style=\\"display:none\\"># S+ i2 ]$ j! O) j  p</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">* i% }; ~3 ^/ h) n</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\\">6 }5 {2 v6 ^3 C. l</font><br /> };<br /> <span style=\\"display:none\\"># W2 |0 e6 \) P+ r</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641
3 Y" g8 ?# w# b" k- |5 F# M6 M
& ^: a9 d- p4 bForcal比Matlab稍慢些。
) H3 {( x7 F5 O+ }& p3 A  F2 L: r, k3 p
－－－－－－－－－－
2 a% F$ D1 Q# U/ W- A8 T; y
/ _+ W8 _/ l: j  j* @- L! H- `再看循环效率。

! \5 _1 t# P( T+ Q" m4 E! RMatlab代码：

1. tic4 c0 U' K3 B* D* h9 H
   
2. s=0;
   2 T) H8 t4 u) }$ p  w
3. x = 0;: L% ~3 I' X% P
   
4. y = 1;; J8 v' g! \: J* r. {2 {
   
5. while x<1    , y4 B) o* F( E+ i6 t) o' T
   
6.     while y<2;        ) R) A- t# H4 ^: a
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));- P\" i7 |3 r) U( \! ^( ?( b
   
8.         y = y+0.0009;        
   : b\" Y1 [- ?3 W( _, N  F
9.     end
   0 ^4 k+ V( L& W% U% s
10.     x = x+0.0011;
    # K/ h# z6 ]0 r) o1 d4 B4 I2 I
11.     y = 1;\" P4 ]. r/ e5 e+ l/ K3 @
    
12. end
    . f/ w5 I7 P. ?( X! V* E+ N! N
13. s. b0 |% B  J, n0 Z( d' P
    
14. toc3 y- i3 s4 f$ r8 c, a* a, U
    
15. ' b0 _7 J) [. ^: A: D& r, A1 O
    
16. s =  E% x- X8 b0 w& t
    
17. 4 F: c# V2 ]+ W- ?
    
18.   1.0086e+006
    / W\" `8 V4 p3 V) c
19. & K0 S6 Y2 ^$ ^: K
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   & K& L% F: Q' S! ]: T* I* y+ [6 ^
2. t=sys::clock();! {# u$ ^! I% t  p( m
   
3. s=0,x=0, ) Q1 h/ q7 A' W
   
4. while{x<=1, //while循环算法； ' U  h' s% y% r- T
   
5.     y=1, 6 u, Q8 S5 O& V\" M( p2 L
   
6.     while{y<=2, & c) H7 C; V8 v9 q2 D
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), % E: Q; b  S) P# [
   
8.         y=y+0.0009
   ) @+ G1 ^) _$ |
9.     }, \" A+ n+ f  _& C  U5 P. L$ m% N
   
10.     x=x+0.0011 3 y; N  s. C6 r8 m4 m  p: V
    
11. },
    + A) g* j2 k4 M
12. s;
    : ]2 T; C$ n0 X/ Z5 V\" Q3 C( e7 s
13. [sys::clock()-t]/1000;

复制代码

结果：
/ z- _- S9 W5 h& H, h% G1008606.64947441
0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

; D6 R* G7 K2 c; I# z1 E6 _－－－－－－－

' e# P. w% x$ H: b% V& lForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   7 T1 i% x5 V  _) Z1 u6 b& r# Q
2. t=sys::clock();
   4 O5 K: M4 P4 ~  N7 |# |; }
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" }: o: N$ i* O
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   1 a6 v1 D/ _% N\" M  x
5. [sys::clock()-t]/1000;

复制代码

结果：
2 i3 O; [* Z. L; O( K( M/ Z: q1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
; w  V1 L) ]$ ?* N, w* y
- P+ |7 o9 U+ G. M( h; j. F/ Zmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 [) B! d* F% H5 e( l$ x0 o
2. tic, |( P0 n/ `- k& {
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);/ q+ u4 X# g2 g5 _1 D
   
4. sum(sum(arrayfun(f,x,y)))
   - V8 C2 B2 Y5 [7 M: [. W3 N
5. toc
   8 p+ ]0 i6 K- i\" m& T  D
6. 2 e, F- ~8 F8 I- y, ]\" |$ V
   
7. ans =  [: Y0 P9 y' E/ @% j
   
8. 7 y\" H0 e7 `0 k: V; Q5 L
   
9.   1.0086e+0069 G+ M! \  [8 M+ _0 F8 w6 {
   
10. 8 C) b. E. N* G( E
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];1 k; u0 ]+ t: E) d- j
   
4. mvar:, h4 l4 C\\" s! Z/ [\\" q2 T
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));   n' f\\" u+ Y$ m5 b' g
   
8. t=clock(),( P# U; J3 T2 k& Q% a
   
10. oo{
    
11. - v0 }* R6 k2 z: p8 g4 W2 V
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],& e( P! F- P( q  w
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 9 {/ D6 Z! O) P4 M; y+ }0 y1 a! R
16. };& c6 j\\" u7 K2 T1 Q
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">1 k; u0 ]+ t: E) d- j</font><br /> mvar:<font class=\\"jammer\\">, h4 l4 C\\" s! Z/ [\\" q2 T</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">  n' f\\" u+ Y$ m5 b' g</font><br /> t=clock(),<font class=\\"jammer\\">( P# U; J3 T2 k& Q% a</font><br /> oo{<br /> <span style=\\"display:none\\">- v0 }* R6 k2 z: p8 g4 W2 V</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">& e( P! F- P( q  w</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">9 {/ D6 Z! O) P4 M; y+ }0 y1 a! R</span> };<font class=\\"jammer\\">& c6 j\\" u7 K2 T1 Q</font><br /> [clock()-t]/1000;

结果：
8 V; h% f6 {; c! A, Q. F7 ]* A  `1008606.64947441
% R; G0 H$ g3 u# N9 S& t% l0.735  秒
3 X, @' X7 H) q# U
7 g) U8 x7 B2 N3 l9 n可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
, C8 R0 K! @7 d8 z
- A' Z& l  ~/ n, [: E2 d; U: {从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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