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

forcal

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

1. //用C++代码描述为：. k: P+ P4 c; d: ^' X6 `
   
2. s=0.0; & N\" {7 m* M- q9 W- u5 b+ m* |9 p: _
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   ! Q\" J$ q* ^/ p
4. {
   8 s, @  ]% j5 R$ z\" K7 I! w0 Q
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   ) _0 e4 Y  O' I' h
6.    {
   % v0 N- o+ [5 I% W' K% u1 Z1 a
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));0 w& U) K& ?  u/ S' X9 E8 Q0 b. I
   
8.    }. S' D4 V3 _\" K4 F
   
9. }

复制代码

Matlab代码：

1. tic
   & s' C8 m; s; k7 J5 \
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);. b3 D2 K, J/ L/ {! G
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))2 {* c5 \1 p! x- n( i2 F
   
4. toc* y2 }& y. D% M: \
   
5. 
   8 b1 L6 G) _7 s2 b; ~\" A# C
6. s =
   $ I; l\" P7 B' Z5 J* h+ }
7. 
   * Q( [; C* Y: g1 e$ d
8.   1.0086e+006' e, d& I! @) y+ a4 N2 I5 U
   
9. - p7 k3 _\" X6 j' v8 q8 t
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];  H# F' E\\" Y7 \7 }! g/ h( @
   
4. mvar:6 p- L% r' m( ^5 b7 h% U. e: C: E
   
6. t=clock(),. R' V6 q6 Q# d/ {$ H
   
8. oo{$ ?/ A8 ^6 g$ _3 J
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 1 g* p) l* }' @& 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]' J  m; E3 |: A+ E
    
14. };# j- e6 }( O% D2 b- Y
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">  H# F' E\\" Y7 \7 }! g/ h( @</font><br /> mvar:<font class=\\"jammer\\">6 p- L% r' m( ^5 b7 h% U. e: C: E</font><br /> t=clock(),<font class=\\"jammer\\">. R' V6 q6 Q# d/ {$ H</font><br /> oo{<font class=\\"jammer\\">$ ?/ A8 ^6 g$ _3 J</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">1 g* p) l* }' @& 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\\">' J  m; E3 |: A+ E</font><br /> };<font class=\\"jammer\\"># j- e6 }( O% D2 b- Y</font><br /> [clock()-t]/1000;

结果：
9 J3 p% ?7 \4 `5 \- y/ Q1008606.64947441
0.641

2 Y! ]0 ~9 ?" ~+ j9 \Forcal比Matlab稍慢些。

/ z, G* x$ b2 r－－－－－－－－－－
. ~! Y4 }3 G6 a& e$ G% t4 T
# Q; \4 h1 |* j1 C再看循环效率。
7 _0 G7 a4 a" x% P: e/ b
Matlab代码：

1. tic: _0 C. Z  N7 ]  T+ i5 z# C. m
   
2. s=0;
   # t5 B- R8 m  K\" f/ a! A
3. x = 0;( c! a# E6 H9 t
   
4. y = 1;& v. u; `4 u% s( j7 }! Q\" `
   
5. while x<1    ; L1 L/ _2 W8 e/ s; }/ Y$ r
   
6.     while y<2;        8 [' X$ A) \) S0 w. i
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));: ]5 W5 L1 q' J2 X* A
   
8.         y = y+0.0009;        : n) E& V7 v2 S
   
9.     end
   & |5 p' ?' g6 C% s$ P8 r
10.     x = x+0.0011;) d! h; h4 w: c- p' X: f( b
    
11.     y = 1;( n; q! ]1 s$ b4 s, X9 H3 D
    
12. end
    2 N1 S\" C) t7 I5 ]1 g& P4 D% ?1 l
13. s
      z3 N9 E/ Q) d1 o/ {2 ~\" K3 M7 R' o
14. toc
    & v( Q7 Z- s# Y, t
15. * o\" q0 }' N' i\" Y
    
16. s =1 }$ [: j1 y! P9 K5 b  Y/ \# ^
    
17. + u% m' y4 A: d6 q: I/ s: P\" @
    
18.   1.0086e+006
    7 `$ X1 k0 f$ O8 g
19. \" I6 a' h5 d8 T) A1 y& P\" x+ P
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   9 E' V; a$ m8 [9 c3 @9 s# ?7 w
2. t=sys::clock();) G( X\" d9 `# v: J! Q
   
3. s=0,x=0,
   1 m0 Q( I8 @+ U0 ~
4. while{x<=1, //while循环算法； ! {1 D) y% \% V4 x
   
5.     y=1,
   ' F1 Q: D7 L; J2 W
6.     while{y<=2, \" N9 r5 O# \+ }9 e3 s
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ' k. D7 I5 E\" }# L
   
8.         y=y+0.0009
   5 ~) p9 c0 c# l* i! }$ [: ]
9.     },
   8 B4 s/ x# w# `% p/ d. x. _0 z
10.     x=x+0.0011 & \5 ~4 M6 X! @' E1 i0 L
    
11. },
    6 B5 U( R+ }2 j& S0 u0 ~
12. s;; y  x- }. B' L1 n. H: ]  D- ]
    
13. [sys::clock()-t]/1000;

复制代码

结果：
, d1 P7 C  t: v7 r1008606.64947441
& |+ F4 Y+ J, i3 c$ ]' P. F0.734   //时间，秒
- [" ^- u: `5 g
( U6 ^* V$ r8 ^! T" x: [: H我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
1 O; U! b& F5 [2 E
' g, Z# {0 }3 j/ U7 F* E' \Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   6 w% \. b. e1 g% i1 Q7 U
2. t=sys::clock();0 o$ r; b0 S; z8 B
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   2 b, [) c- h2 i5 o
4. sum["f",0,1,0.0011 : 1,2,0.0009];1 Z$ c/ G6 K# A/ D
   
5. [sys::clock()-t]/1000;

复制代码

结果：
) ?- O; K  ^0 x: ^1008606.64947441
0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
) Y: {0 a  @- [1 N. r
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . L2 f8 e4 L* r\" E' J$ W) f7 R  X9 d
2. tic
   \" T! D: n1 P# I# F% P( f) c
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   9 m$ ^4 ]! |. u! s9 v# g( X
4. sum(sum(arrayfun(f,x,y)))6 y6 I: q; D- s/ V8 M
   
5. toc
   + T) D: F4 W# M# ]
6. 
   \" J2 n. d- _& c9 p8 w, U6 Q9 [
7. ans =) v& r\" O* F# M( E0 I5 V
   
8. 
   ! G  c7 @% j+ l6 ~5 t\" r2 {\" Q! E
9.   1.0086e+006
   4 R! y6 y+ d' v+ W2 z
10. 
    3 {, d# ]1 M5 k  R\" T& O% L5 q
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. : O+ n4 e' R# r, @8 G  ?# Z
4. mvar:
   
5. 6 @4 E& \% D9 L+ }& C' K4 \0 t9 |
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 8 B8 c  U9 N/ u7 S3 Z& B. v
   
8. t=clock(),6 d2 @5 x! y6 }6 a- V$ X6 `, i& c
   
10. oo{
    
11. * e8 A6 M# X1 |7 ?/ @: p
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],( }' Y8 I- C5 R) ]6 y) k  U% \1 j
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. ; f6 T; `8 P. i
16. };# K' M& r9 U\\" ?/ E/ o& d+ i% x
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">: O+ n4 e' R# r, @8 G  ?# Z</span> mvar:<br /> <span style=\\"display:none\\">6 @4 E& \% D9 L+ }& C' K4 \0 t9 |</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">8 B8 c  U9 N/ u7 S3 Z& B. v</font><br /> t=clock(),<font class=\\"jammer\\">6 d2 @5 x! y6 }6 a- V$ X6 `, i& c</font><br /> oo{<br /> <span style=\\"display:none\\">* e8 A6 M# X1 |7 ?/ @: p</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">( }' Y8 I- C5 R) ]6 y) k  U% \1 j</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">; f6 T; `8 P. i</span> };<font class=\\"jammer\\"># K' M& r9 U\\" ?/ E/ o& d+ i% x</font><br /> [clock()-t]/1000;

结果：
7 U/ Y' R3 a; X1008606.64947441
0.735  秒
6 }. C" F8 @7 B, W
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

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