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

forcal

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

1. //用C++代码描述为：0 L8 e9 q8 e; ~( N+ h
   
2. s=0.0;
   : R$ t- u+ b) A! q: W* Z
3. for(x=0.0;x<=1.0;x=x+0.0011) $ Z2 r; v8 ^2 p
   
4. {
   7 x$ N! a( c' M3 E! N- `
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   8 D9 C\" e3 |: a
6.    {0 \  B- M' r/ O8 S* H% V2 W; Q
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % \\" p, I4 j; L0 c. m8 K# i
8.    }
   \" {  ]5 ~' W- @
9. }

复制代码

Matlab代码：

1. tic. o* w  ~8 Z; }- w5 Y, x% E: X
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);7 r9 i* s  x8 Z7 L
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   . X$ ^. T5 i0 v7 G
4. toc* M8 _/ z- T2 D! n3 V
   
5. . j0 W& `* l' U% b: Z
   
6. s =& o2 L3 G. {) ^0 r
   
7. , _0 y- w' ^0 \& a& r$ r
   
8.   1.0086e+006\" f- o0 ]) R$ J) m\" x
   
9. 
   0 L/ l# h7 h: l\" a% y1 x
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ) |8 a! r6 M/ ^
4. mvar:
   
5. / B9 W, Q/ t- [, W- G
6. t=clock(),% a% U% g( j6 Q7 ~0 O* e' N! k3 R; i
   
8. oo{9 j% R% W) q3 k& R/ Y' n& u. |
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],9 H, V1 Z' L# O\\" o5 @# X/ I9 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. 3 J: z5 g\\" M1 N6 ]; r- n+ x
14. };
    
15. , [  o5 q0 t% W
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">) |8 a! r6 M/ ^</span> mvar:<br /> <span style=\\"display:none\\">/ B9 W, Q/ t- [, W- G</span> t=clock(),<font class=\\"jammer\\">% a% U% g( j6 Q7 ~0 O* e' N! k3 R; i</font><br /> oo{<font class=\\"jammer\\">9 j% R% W) q3 k& R/ Y' n& u. |</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">9 H, V1 Z' L# O\\" o5 @# X/ I9 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\\">3 J: z5 g\\" M1 N6 ]; r- n+ x</span> };<br /> <span style=\\"display:none\\">, [  o5 q0 t% W</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641

* J" ~! f/ }; E4 ^8 @8 w$ uForcal比Matlab稍慢些。
( `" l& K; v4 R8 g3 {3 |
5 L2 u& C: Q4 @) N1 ~2 q  P+ B! d0 C－－－－－－－－－－
; V% k( i4 d4 Z9 I# f" y1 i' y1 q
2 `. S% f4 i/ Y再看循环效率。

Matlab代码：

1. tic
   ' ?6 u4 N( S: Q6 O5 N5 R: v
2. s=0;
   3 Z1 [7 M) p7 h
3. x = 0;+ v2 s7 u2 N4 K5 n- ~( U
   
4. y = 1;
   2 D# E% @! p/ M
5. while x<1   
   ! [1 ^4 ]0 h) {2 x- x
6.     while y<2;        
   ! E& ?) X. ]8 M9 S5 B\" F\" P+ k( E
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( s$ I& G: E6 u  g, A1 i, N  @
8.         y = y+0.0009;        
   2 l8 `3 r+ o8 m
9.     end: a5 m0 x! [: L& p% ?
   
10.     x = x+0.0011;
    # \9 g; H2 _7 h# a& b
11.     y = 1;
    + Y; K/ x% S! J9 X
12. end
    : p# v) J# R$ q
13. s
    5 ~/ Q0 [% q' T, _6 V# C* c
14. toc8 j' y8 Q  b0 v, E. c  q; l1 O* }
    
15. 3 D' J' ~4 j1 _( F& z  @
    
16. s =
    1 J+ o/ i+ F  ]$ [/ Q; n( H
17. 
    3 s' k) P- ~, Y
18.   1.0086e+006
    % h+ E% {6 b- Y
19. 6 }9 X. E4 A# K% @# l
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   / ?3 Y0 U% a. n. p: g+ N# r
2. t=sys::clock();
   % u' c. d3 }  T2 `! U6 E
3. s=0,x=0, \" J( j6 U( K7 W+ ~4 ~' [
   
4. while{x<=1, //while循环算法； 8 U2 W+ m\" l/ P+ b0 `3 E0 U\" p6 j
   
5.     y=1, ( B1 {* A; n7 y' {& a: i
   
6.     while{y<=2, / `- M  i: x* Q$ h
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   . n2 b( a! E6 b# {6 f* c
8.         y=y+0.0009
   ( X7 Y  w$ N$ Z( ?  c
9.     }, # l6 Q+ X/ c) H$ p\" Z( E! ]7 N; ]
   
10.     x=x+0.0011 ! K- M, O0 R. h
    
11. },
    ( r3 I0 S& L8 J) y
12. s;
    - Z$ K5 i/ L) X% c2 j( }
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
: W7 S2 H, z* W
3 ]# @; v: {4 Q% b+ u/ g% n( Y我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
7 p2 U; v' }9 b$ S3 p
3 z7 a* v, G0 T9 W3 l－－－－－－－

' `$ ^5 V" D: {Forcal中还有一个函数sum专门进行这种计算：

1. mvar:: Y$ a, X/ D2 k' g7 K
   
2. t=sys::clock();
   - X0 d& Z8 z( l* w: k1 }/ ^
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 t! I1 }1 b% x
4. sum["f",0,1,0.0011 : 1,2,0.0009];
     g3 j% E: f4 N
5. [sys::clock()-t]/1000;

复制代码

结果：
0 ?$ q9 R. m& k5 O$ Q+ L1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
" J7 \( r, `& Q/ U0 U; P
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! I2 |( A  G. b/ U  M. A/ O
   
2. tic5 ^  ^  p+ U7 g3 T, u
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   9 q. i- {' V5 Z: A
4. sum(sum(arrayfun(f,x,y)))
   ' H1 |6 P) j0 A1 `$ F. ]& g
5. toc
   $ h- W7 f+ |2 F- M2 i: [, k9 {; `
6. 
   % E0 d6 K& U4 |& f) ]
7. ans =
   4 Y: g0 u* A, s6 t' v0 _4 R
8. 8 A  D+ w8 K- T# P
   
9.   1.0086e+006
   & l' x+ a+ l7 A0 h
10. 
    0 c6 w+ P) G* K% U, Q8 T! L7 E
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. - |% @  g- M8 s4 @! h3 G& M8 ]: @
4. mvar:
   
5. & W* C2 B! X% e! z( L
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));   i# p7 c5 }& I
   
8. t=clock(),
   
9. 2 O' P$ u/ n\\" X: ^/ Y( V
10. oo{- F  F  C+ Z  v
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],+ g- H) g% q  P' S( R
    
14.     arrayfun[HFor("f"),x,y].Sum[]( h4 r- N/ T3 F+ @. O9 ?5 E% O
    
16. };
    
17. ; R( d3 W; `/ [3 G! y
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">- |% @  g- M8 s4 @! h3 G& M8 ]: @</span> mvar:<br /> <span style=\\"display:none\\">& W* C2 B! X% e! z( L</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">  i# p7 c5 }& I</font><br /> t=clock(),<br /> <span style=\\"display:none\\">2 O' P$ u/ n\\" X: ^/ Y( V</span> oo{<font class=\\"jammer\\">- F  F  C+ Z  v</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">+ g- H) g% q  P' S( R</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">( h4 r- N/ T3 F+ @. O9 ?5 E% O</font><br /> };<br /> <span style=\\"display:none\\">; R( d3 W; `/ [3 G! y</span> [clock()-t]/1000;

结果：
" Z- n, t2 T  D1008606.64947441
2 b; Y$ a( b4 U3 H0.735  秒

5 J0 ^- F  ~, T- Y. c1 K可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

2 q0 [% m5 @' a' S3 d+ A" r－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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