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

forcal

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

1. //用C++代码描述为：
   : r- `+ L; [) D: |/ j
2. s=0.0; 8 u2 I! K( C\" F
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   8 q, o2 K) ?$ H. Z- q
4. {/ W3 L1 o7 @9 q4 Q
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)/ H) f1 q1 M2 a7 k  A& B5 s
   
6.    {
   , m- s% x2 a  c% z: e6 h
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));0 r& N# b4 Y\" h* c! S$ N  f
   
8.    }* T0 b! s: e( g* U. ]
   
9. }

复制代码

Matlab代码：

1. tic: G% A+ W' }! `# t% C: j
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   2 W0 w1 U$ m7 C; p* {\" \, B! ?\" N7 ]
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   0 L: j+ `& w( _$ z. a4 A4 K% }0 i
4. toc. L$ T; q9 C' _) m  }) S8 e
   
5. 9 M% u4 U/ C$ u* \% g! D
   
6. s =
   * j# \/ w1 |) V# l8 Y
7. 
   6 M( I' n/ {) C( \' t, H
8.   1.0086e+0063 g. a  p, [& L+ z# O
   
9. \" }% K$ R\" a! L
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];0 H  Z5 Q/ l) w$ k- J
   
4. mvar:
   
5. % h- \* z  ^, u: }4 B& g7 l
6. t=clock(),, r+ M\\" r2 A6 U5 M
   
8. oo{5 h1 @  D: A. R
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],4 b# y* Z: }2 s* J: {
    
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. 7 m- m+ m( L\\" |! U
14. };) R\\" M\\" ~3 v- W! O, R6 h, k5 J3 Y
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">0 H  Z5 Q/ l) w$ k- J</font><br /> mvar:<br /> <span style=\\"display:none\\">% h- \* z  ^, u: }4 B& g7 l</span> t=clock(),<font class=\\"jammer\\">, r+ M\\" r2 A6 U5 M</font><br /> oo{<font class=\\"jammer\\">5 h1 @  D: A. R</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">4 b# y* Z: }2 s* J: {</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\\">7 m- m+ m( L\\" |! U</span> };<font class=\\"jammer\\">) R\\" M\\" ~3 v- W! O, R6 h, k5 J3 Y</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.641
7 q* T. ~2 w9 ~3 o0 L& C5 }
1 a: `, K  N$ v& O0 J+ fForcal比Matlab稍慢些。
/ w  G/ b: L8 E6 c, W
/ F  Y% p8 O7 Z( P6 @' N4 v－－－－－－－－－－

( _$ E. I' C5 \% N7 D. I再看循环效率。

; j- [: G; _4 \: s% i& ZMatlab代码：

1. tic  ^  c: D\" c. b9 K
   
2. s=0;
   3 O! N, v1 m\" ?5 I! l9 X3 R4 Q
3. x = 0;/ f' z1 P. I7 t( F( h+ _  {
   
4. y = 1;7 M4 L/ E) Y5 w; m
   
5. while x<1   
   6 ?4 ?. M4 C9 Y\" g& n: t
6.     while y<2;        : g! e, j) V- p
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));7 o# n7 J) J* }\" c\" a, ]
   
8.         y = y+0.0009;        
   # i! B3 ^- }1 B5 y8 o) s
9.     end
   ' r( d0 v' x6 r) K  x! x) T
10.     x = x+0.0011;7 S' Y3 c) u5 L; r\" t
    
11.     y = 1;
    8 L0 X+ Z, h8 O' e
12. end
    $ r2 k9 q5 y% O! k+ |( |
13. s
    8 u) ]- ~# p+ R, m
14. toc$ j4 v! ^( @) @# U  [, X! F
    
15. 
    \" ~' I7 k# f* y3 l- u) {
16. s =
    % }; p; w+ @4 m' Z4 P, W% q
17. 9 y' ~& H7 c# b2 _$ g
    
18.   1.0086e+0062 f; f\" h8 a8 I( I5 e0 H; P
    
19. 9 }9 j; h2 |4 h( Y9 P( O
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   1 v  Q, Z  ]% O; v; ]# }* E
2. t=sys::clock();
   5 V3 P) F/ I( N0 {; t
3. s=0,x=0, 1 Q\" s  v! M; j4 a' z
   
4. while{x<=1, //while循环算法； 0 |3 E7 q( Z  s9 o
   
5.     y=1,
   . A; v/ \+ G% d% p$ Z# f* F
6.     while{y<=2, 9 A: V: `' h' c) d. h
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), & F: C- g: M1 E
   
8.         y=y+0.0009 ! e( a# w! `\" Y1 S
   
9.     },
   + ]1 [1 M' B) G
10.     x=x+0.0011 9 B6 I\" H5 w; ?
    
11. }, . d- m' S: Z/ E! L
    
12. s;
    6 ?- E9 Z1 r7 r4 G
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

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

3 x! B, _& w3 N( m) [1 Q－－－－－－－
2 {/ p5 O9 h, x! X
2 t" X" B& w2 Y  N( wForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   $ J1 c9 i/ k4 T2 J7 L
2. t=sys::clock();! g: V( P, |; b; S
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 2 P: h9 V0 [1 X) k8 h/ J+ ^$ t; V. }
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];/ `# o+ D\" j5 [0 _5 c
   
5. [sys::clock()-t]/1000;

复制代码

结果：
( ?, N5 w$ _6 g! X  q1008606.64947441
2 ?3 \& F% ]1 s  |1 `6 ]9 y" h0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。

& p- G3 L- }$ F) A8 Tmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));9 }! m3 o2 f: |8 C. h9 e- T( x
   
2. tic, C$ |1 w2 Z# C8 w! n. D
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);3 }4 A5 }  K6 X4 Z1 u
   
4. sum(sum(arrayfun(f,x,y))), {( K$ _& Q/ _2 W\" j; H
   
5. toc, n0 ]/ t* V; A! F
   
6. + Q( j0 k& ?+ R; H0 p9 Y9 ^; W1 w
   
7. ans =
   + H\" g- q4 W5 n! L$ U$ j1 c' p
8. 0 t* h5 k! Q% P5 `5 d
   
9.   1.0086e+006( `. e; s# S+ v/ X
   
10. 
    8 ?! c% x6 }; ^+ u3 t
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];0 X\\" o5 [' I9 L3 V! M0 T# o7 m
   
4. mvar:4 a% }; Z3 b9 W: Y- I4 w& d
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));   ^- N6 E4 z\\" c
   
8. t=clock(),3 \, r! k1 M7 D) O3 I) Y6 Z
   
10. oo{, U- M! X8 Q. P5 p
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 5 t9 V+ Z3 F* ~
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. 7 K9 J; [9 T3 n5 s) u- o! J( @) r
16. };7 e4 u1 B\\" H+ ^5 ^
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">0 X\\" o5 [' I9 L3 V! M0 T# o7 m</font><br /> mvar:<font class=\\"jammer\\">4 a% }; Z3 b9 W: Y- I4 w& d</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">  ^- N6 E4 z\\" c</font><br /> t=clock(),<font class=\\"jammer\\">3 \, r! k1 M7 D) O3 I) Y6 Z</font><br /> oo{<font class=\\"jammer\\">, U- M! X8 Q. P5 p</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">5 t9 V+ Z3 F* ~</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">7 K9 J; [9 T3 n5 s) u- o! J( @) r</span> };<font class=\\"jammer\\">7 e4 u1 B\\" H+ ^5 ^</font><br /> [clock()-t]/1000;

结果：
! i2 L% }9 t4 L: m7 ?& g/ Y) F1008606.64947441
0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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