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

forcal

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

1. //用C++代码描述为：
   ; g9 a$ W2 V2 ?9 v
2. s=0.0; 7 ?3 B; U) \1 k4 `6 L
   
3. for(x=0.0;x<=1.0;x=x+0.0011) 8 a0 C9 P7 A( c: t* l
   
4. {& j- h* |9 z* H7 g( l# c/ i
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)* J+ u9 F) ]& c1 |
   
6.    {
   8 q; q/ ?4 Z\" j- C  i' s' I
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 [& O9 x' p# e. q; R
8.    }8 n3 w# n; }- T. S/ e& G
   
9. }

复制代码

Matlab代码：

1. tic/ q. J\" d& f. H5 k1 I' P' D
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   & _3 H2 R4 E6 b
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ' ~* {: @. R2 K  H
4. toc
   & l: W) n0 K' d0 o) a; @+ ]
5. 
   5 C# \3 j* P% F  L
6. s =
   - j* P- T0 w7 @* b' G% R/ k) |6 u
7. 
   5 S8 C: a\" p: X- o
8.   1.0086e+006% Z, I4 V1 L) v5 E! r
   
9. 2 ~+ L( H& G; A
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];3 ~$ B: \4 y/ C$ k
   
4. mvar:- z# ~% d) o  N' B* o& \$ {1 Q
   
6. t=clock(),
   
7. 8 T3 ^: A: U2 ~  e2 x& h
8. oo{5 X& F0 R) Y7 r; C1 K
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. ) j# J, @( @3 h( B; ]1 N, p
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]: H+ Q# G. x+ ^
    
14. };
    
15. + N5 q* {3 r% o5 Q7 a; P
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">3 ~$ B: \4 y/ C$ k</font><br /> mvar:<font class=\\"jammer\\">- z# ~% d) o  N' B* o& \$ {1 Q</font><br /> t=clock(),<br /> <span style=\\"display:none\\">8 T3 ^: A: U2 ~  e2 x& h</span> oo{<font class=\\"jammer\\">5 X& F0 R) Y7 r; C1 K</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">) j# J, @( @3 h( B; ]1 N, p</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\\">: H+ Q# G. x+ ^</font><br /> };<br /> <span style=\\"display:none\\">+ N5 q* {3 r% o5 Q7 a; P</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641
8 i5 ^  w. L$ c9 U6 K
Forcal比Matlab稍慢些。
# B/ C+ R5 ~+ D- h
6 \( q4 b1 y5 Q" }; y4 ]; i$ j－－－－－－－－－－

+ p5 T- D) z5 a+ b+ s2 o/ W再看循环效率。
5 D9 P7 @  R8 q" ]
0 L1 [, R( C  F5 }Matlab代码：

1. tic( c) H7 ]( N1 \' x( ]. j
   
2. s=0;
   % h. O0 E& }) i- U7 n+ j! l7 i
3. x = 0;2 `! B& s\" `# x2 S, O5 |- E+ z
   
4. y = 1;* m' O, d9 P& s  c# d! d
   
5. while x<1   
   0 A  Q. y0 ~5 C) J/ H4 C
6.     while y<2;        0 b4 \\" e: N) k6 \
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; `$ v* x3 e; x5 `5 A; F5 f1 l\" x0 f
8.         y = y+0.0009;        
   ! T2 j; N; ]! M\" |$ J. v
9.     end) {% S0 R! g2 x8 r5 g
   
10.     x = x+0.0011;
    & K9 i\" z; t& V* ?! ?
11.     y = 1;
    \" H+ f  M5 v1 b7 g+ [8 F
12. end5 `1 v8 |, p% O* T: V; c
    
13. s
    $ K$ Y2 w1 H9 |2 d( R
14. toc
    / t0 |2 }) A/ a! ]8 ~/ i8 N
15. 2 d1 v9 Q( _5 o: `1 Q  Z
    
16. s =
    ; |\" y) b3 B- ?0 J9 P( r( v
17. $ Q% S. h/ a2 n4 h9 A9 x2 c1 `
    
18.   1.0086e+006\" Y3 s* s! N. V/ t
    
19. 7 X- ^. S\" s! g
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:( _  l9 [5 a# p. Q) e  k
   
2. t=sys::clock();
   & V  M- e# R) K$ y5 A% q
3. s=0,x=0, ( W* _. f. \# V
   
4. while{x<=1, //while循环算法；
   \" h2 N& E3 V8 m\" ^# u
5.     y=1,
   - B' ]) S2 L\" g$ z. q9 \5 o
6.     while{y<=2,
   , J$ H+ o) X  h% A\" Q0 M4 B
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ) x  Y( Y! d2 Q' m
   
8.         y=y+0.0009
     F& Y+ B9 v\" h& Y) j
9.     },
   5 c+ ^1 M& E5 q/ i! K3 Q2 }
10.     x=x+0.0011 7 U, J: u2 F4 q0 T  l) f& m$ _
    
11. },
    : r* t3 P4 L  f
12. s;' R( n' |! C2 `; n! y
    
13. [sys::clock()-t]/1000;

复制代码

结果：
, ~& H4 {9 {5 B/ M5 N1008606.64947441
, g8 S( N( U' R/ }9 Y0.734   //时间，秒

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

－－－－－－－

" J! j5 u! D3 T' F' M; a$ OForcal中还有一个函数sum专门进行这种计算：

1. mvar:+ H, Q6 W9 ~( B. P& \1 U0 h5 W
   
2. t=sys::clock();
   0 a; f- f( v7 z\" d# k
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   $ w8 l* H9 c* n6 W
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   5 K& ?+ D: p) Z3 Q* t
5. [sys::clock()-t]/1000;

复制代码

结果：
' g0 Y6 x' q, c! _1008606.64947441
7 ^* E- y: `5 ]  L  t0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
: h! B" b9 M  A2 E, W) v
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( R( B( ~  ]9 J& d4 ^/ [) X1 }
2. tic2 A6 x- d* H( Z& }) P8 ?: B) k  q
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);1 A( V; o0 \2 l# t5 \6 [+ Z) `* Y8 ~3 j
   
4. sum(sum(arrayfun(f,x,y)))
   * y\" g1 p& ?9 @: \  Z% u1 ]\" K
5. toc2 L2 `( T4 A1 g: ~/ _7 J. _3 q
   
6. 
     d5 }# L7 F! l
7. ans =3 Q0 V# }7 C) ~  |, K4 G0 N
   
8. : t. a: k. [4 e2 F
   
9.   1.0086e+006# @+ {* D$ x# V( t\" {
   
10.   h/ f- h$ \# M! g8 ?
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; W( }' _( x/ R9 e8 N9 g
4. mvar:
   
5. 3 z; G  @* \6 Y  Q
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 6 g1 u8 g3 C3 N
8. t=clock(),
   
9. $ l: d; w/ Y- x+ F* s* v1 ]
10. oo{
    
11. 1 K  @! t+ t1 f, }6 l0 S, ~* Z
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( V: [* W& K5 K' @
    
14.     arrayfun[HFor("f"),x,y].Sum[]# E' _* {  u, R' R  ~, W
    
16. };
    
17. % r: B5 k, o' ?+ C8 [/ Y
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; W( }' _( x/ R9 e8 N9 g</span> mvar:<br /> <span style=\\"display:none\\">3 z; G  @* \6 Y  Q</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\\">6 g1 u8 g3 C3 N</span> t=clock(),<br /> <span style=\\"display:none\\">$ l: d; w/ Y- x+ F* s* v1 ]</span> oo{<br /> <span style=\\"display:none\\">1 K  @! t+ t1 f, }6 l0 S, ~* Z</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( V: [* W& K5 K' @</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\"># E' _* {  u, R' R  ~, W</font><br /> };<br /> <span style=\\"display:none\\">% r: B5 k, o' ?+ C8 [/ Y</span> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
( v. `& \+ y; w: i( J
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

1 y% F0 d1 B- Z: X3 x& R% j- U－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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