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

forcal

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

1. //用C++代码描述为：3 L6 u& F, L6 `9 o$ s. N
   
2. s=0.0;
   8 A; T2 G% L  B! n4 C, g: k
3. for(x=0.0;x<=1.0;x=x+0.0011)
   5 v3 x0 b; }( a4 o& x# e
4. {$ S$ f# t4 x7 ^
   
5.    for(y=1.0;y<=2.0;y=y+0.0009), o1 Q- ~- N' v* t2 B. y5 M9 N
   
6.    {6 B7 O& A& y, O. f( i2 K. a9 ~
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   : w$ b9 o& [' ~6 ~2 ^9 w
8.    }\" s' _' e, ?' Z
   
9. }

复制代码

Matlab代码：

1. tic# }* C; p  I3 ]& _! ?
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   1 c( h/ k2 `\" c& L5 q# A
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   ' ?1 L4 i, B5 J# ^& ^
4. toc
   & ?: r$ m+ D/ b! l% n; f9 c4 a' `
5. ; Q6 c& ~2 j( B7 `+ I* S, k
   
6. s =
   5 N' G5 X4 n2 b+ g4 Q9 ?& T
7. 
   \" X$ W4 V5 L( S/ p$ T7 ?
8.   1.0086e+0066 c! p6 g\" l# `
   
9. 3 H' i1 D% Y% L) M1 C( o
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 2 g$ V# w: w- n  |3 }7 D/ S( Z/ M
4. mvar:
   
5. / F7 ?# T5 s$ H  m! F- s
6. t=clock(),
   
7. & w, `. S5 Q) e# d
8. oo{
   
9. / ~7 M2 T8 |/ H5 V3 x
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 8 C6 |' ]  q+ ]5 W% |
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. 2 W+ |- L* {9 Y$ x  T0 t2 W
14. };7 b3 r* P$ Q% b' [2 T- [
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">2 g$ V# w: w- n  |3 }7 D/ S( Z/ M</span> mvar:<br /> <span style=\\"display:none\\">/ F7 ?# T5 s$ H  m! F- s</span> t=clock(),<br /> <span style=\\"display:none\\">& w, `. S5 Q) e# d</span> oo{<br /> <span style=\\"display:none\\">/ ~7 M2 T8 |/ H5 V3 x</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">8 C6 |' ]  q+ ]5 W% |</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]<br /> <span style=\\"display:none\\">2 W+ |- L* {9 Y$ x  T0 t2 W</span> };<font class=\\"jammer\\">7 b3 r* P$ Q% b' [2 T- [</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
5 ~/ k6 k; v7 u8 U% R! U  y0.641
2 I- H# y: l  O8 ?$ y; f1 E$ f$ a; `
# v- B# Z( ~0 x1 h. m) K4 i7 K" AForcal比Matlab稍慢些。
; _% S5 W2 J! P; ?/ h# i
8 n& ?: @) ]2 l! G0 {/ r－－－－－－－－－－

/ O* j/ u. Y, Y# }8 M) r再看循环效率。
+ q- k% c3 I2 T
Matlab代码：

1. tic; s( e: n* g6 q
   
2. s=0;
   ) ~5 W( ~9 o7 k# H* W
3. x = 0;' a. Z$ {, K6 k; C
   
4. y = 1;
   2 h' @1 G2 ?/ p4 r- j5 |
5. while x<1   
   7 r8 T  f6 W6 ^- [5 v8 n
6.     while y<2;        
   : o) {% p  Q$ [
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - |% E: v% t9 i
8.         y = y+0.0009;        ! T7 @( u% X* m. F
   
9.     end
   9 p* E- e- P/ h4 Z/ Y2 \! H/ d% R3 d
10.     x = x+0.0011;
    7 m% O8 l# u7 ^& V$ [$ m0 e
11.     y = 1;: I# t8 }' v1 D. ~% z$ a6 g
    
12. end
    ; |- Z5 j8 O. z( U' R- j
13. s0 H- z! n' @2 w- O9 A
    
14. toc
    - `; y3 P* c' Y) @2 T3 A; p8 C
15. 
    7 _2 I+ m- B\" A7 Y' s! [
16. s =
    / Q5 X5 c5 a% H
17. 
    9 d\" Q- T) \! I* M1 I, ?
18.   1.0086e+006
    , f/ J* y: i4 x$ x' c
19. 6 l1 R+ P( G2 k2 s1 x( T9 Q& V
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:3 Z4 l  N. {0 h. b6 P5 r! w, u: c
   
2. t=sys::clock();8 i- d; l7 V9 ^/ d& I# h
   
3. s=0,x=0, - a/ x3 ^' e# o6 }* Z7 p* d& V
   
4. while{x<=1, //while循环算法； 0 Q7 I% Q& Q) P; e( L) [3 J
   
5.     y=1,
   8 U3 y\" T8 u  S\" m* ^1 m0 W& E
6.     while{y<=2,
   1 K9 f( J9 ~- C' l
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 0 a8 N\" d\" \7 |% n- q
   
8.         y=y+0.0009 # H8 Y7 F3 `2 M: t- z: W6 r
   
9.     },
   ( J% E5 V; K+ R1 n
10.     x=x+0.0011 7 P# d; J* v4 o$ w$ b/ v
    
11. }, 6 k: x( b9 U. _  [& U
    
12. s;4 S: `0 M0 a  C6 y% m! Z3 a
    
13. [sys::clock()-t]/1000;

复制代码

结果：
& t, g  Z/ q( a9 E# _- t  D& {2 D9 m- `1008606.64947441
5 e7 K% ]! y2 |" n7 K8 A# ?$ J' e& p0.734   //时间，秒
" a5 q) s& U( H2 d) L- R
/ {/ {. w2 Y* Z: o2 x* a我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

( _, D/ m- i& c0 {- N7 o1 p2 U5 N7 j5 O－－－－－－－

" t! J  s/ p, @# A3 fForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   ; ^3 l# l5 G\" I0 l5 {1 P
2. t=sys::clock();4 F9 \+ D5 C! |, _
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));   Q! ~$ f+ {  @% s6 @4 y
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];5 t4 z% N0 e\" K+ o! Z
   
5. [sys::clock()-t]/1000;

复制代码

结果：
( W; i- K4 \' K9 `* S6 `1008606.64947441
8 o+ o" C% K- K  X# D' @0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
' d7 d! z  x; |# g0 e& n7 M
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; C- n( p7 h) Z) p
2. tic
   7 C* F8 o$ e\" q' r7 J4 R
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   # R4 X, _! x* h! l# p6 z, i; j: P7 l
4. sum(sum(arrayfun(f,x,y)))
   6 a\" F% j! J; d# v- q7 O1 Z2 r
5. toc' p, e2 M7 h! s+ V\" p
   
6. ; W# q/ d) t' q+ L
   
7. ans =  C  p1 _- v7 w
   
8. 
   + j( W+ @0 O; K) D2 A2 F+ W) U0 y
9.   1.0086e+006/ b- [4 _) x' `\" y6 I2 d
   
10. 
    5 a! @- ]6 R, f5 J! A% D2 f' N9 I\" c
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. $ l\\" A0 d\\" ~6 c& J+ W2 m) e, x, ^
4. mvar:
   
5. & V# v6 {/ P+ Q
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); & P/ Q7 x# U4 ^
   
8. t=clock(),3 Y: K\\" n' m) g* t/ n
   
10. oo{! q! Q- K; x' K$ {+ l3 v
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. $ r. H6 Y' T! t  V0 `: @
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 3 D\\" y  a9 ?9 A0 s. |7 w
16. };
    
17. % m9 i/ }+ ^% I
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">$ l\\" A0 d\\" ~6 c& J+ W2 m) e, x, ^</span> mvar:<br /> <span style=\\"display:none\\">& V# v6 {/ P+ Q</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">& P/ Q7 x# U4 ^</font><br /> t=clock(),<font class=\\"jammer\\">3 Y: K\\" n' m) g* t/ n</font><br /> oo{<font class=\\"jammer\\">! q! Q- K; x' K$ {+ l3 v</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">$ r. H6 Y' T! t  V0 `: @</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">3 D\\" y  a9 ?9 A0 s. |7 w</span> };<br /> <span style=\\"display:none\\">% m9 i/ }+ ^% I</span> [clock()-t]/1000;

结果：
1008606.64947441
7 _5 c( L  l8 D& l& y$ v# Z0.735  秒

. L' E: J" |& ]* i可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
! Y- O; }7 S% ~/ s. n
, k7 H+ |: J- m7 D& r0 ]! a从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。