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

forcal

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

1. //用C++代码描述为：* N: P/ i2 m5 t- Y& W  ], R
   
2. s=0.0; - B' R; z\" k) W( K3 B. D+ ?
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   # q. j) @( u; |! j% `
4. {
   / j- n, X. J\" M) [' [
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   : r+ a/ U$ D+ [4 T
6.    {# f( W9 o! a) J\" M0 Z8 C8 U
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   # O, w# m( d\" p6 J' M2 O2 c
8.    }2 y' r; W# i+ H4 I- w
   
9. }

复制代码

Matlab代码：

1. tic
   * F  R# ^( j, T6 P# e9 i
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);! Y5 X9 I& G) ]2 y& W# h
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))), m- _5 T2 |/ [1 y+ D7 y% W
   
4. toc
   \" o6 O# c+ e; v
5. : q7 {1 [6 r& N\" }0 m
   
6. s =
   / R9 w2 Q( S( D% V9 f+ d. m- f
7. 
   * V1 C! F/ v3 Y
8.   1.0086e+0061 F+ u0 _5 P\" J2 |2 }' C8 O, V
   
9. \" k3 q# d6 q$ g9 G7 N
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; m) k' a3 K7 |: C$ o, \! j
4. mvar:
   
5. ; Y2 v3 ]8 E) Q1 L8 Y
6. t=clock(),4 F! d  P\\" X; F8 e, F* J: t
   
8. oo{
   
9. / ]# t# ?4 U$ O\\" `
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 2 u: Q7 C4 L6 [! f6 I7 ]. `
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. : f7 P7 |* X+ n. a9 a
14. };
    
15. 7 _! P9 r- c8 \3 p
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; m) k' a3 K7 |: C$ o, \! j</span> mvar:<br /> <span style=\\"display:none\\">; Y2 v3 ]8 E) Q1 L8 Y</span> t=clock(),<font class=\\"jammer\\">4 F! d  P\\" X; F8 e, F* J: t</font><br /> oo{<br /> <span style=\\"display:none\\">/ ]# t# ?4 U$ O\\" `</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">2 u: Q7 C4 L6 [! f6 I7 ]. `</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\\">: f7 P7 |* X+ n. a9 a</span> };<br /> <span style=\\"display:none\\">7 _! P9 r- c8 \3 p</span> [clock()-t]/1000;

结果：
7 a" }/ D1 w  T; }1008606.64947441
; d: O2 I2 c8 P' Z# ]( e( r3 i0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－

再看循环效率。
4 y1 N! I5 @. p( a' p& d& A
Matlab代码：

1. tic
   - w+ B# c\" ^+ J% K9 M
2. s=0;
   4 i6 F' q3 J7 ~9 A
3. x = 0;7 {$ B- U9 q5 O; B0 y
   
4. y = 1;) \+ d9 p9 I+ `
   
5. while x<1   
   & M9 n9 g) m) u5 [- s0 w& h  N/ L
6.     while y<2;        ) J) W' o  Q# b\" K2 y% [1 C
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! h$ F$ X8 ~\" S, c, K. j( _* \4 {\" H
8.         y = y+0.0009;        
   # I/ y. U' s- a: T3 g% M6 O
9.     end
   - B% C6 F5 f& f
10.     x = x+0.0011;
    7 Q4 Z. Q' T\" ~7 c# N\" {% I
11.     y = 1;
    $ ]! p! [& D# Z\" l- ]4 J
12. end
    % s9 N6 n5 h& ~; I9 `
13. s
    / A/ V+ ?/ O. w& f- S# p\" C5 o( m
14. toc, l# _. o$ D! b. c' U. y
    
15. * A, I# Q% A5 p+ s7 M
    
16. s =
    6 w) g\" d9 I: \\" _
17. , u5 j: p& E0 l& r
    
18.   1.0086e+006% B! @8 T\" M: _+ X
    
19. 4 \8 L& Q( W; z& O1 J5 `
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
     ?; A; v/ u* Y3 |# e6 I
2. t=sys::clock();4 _% y& W( Q( _; F& n7 R, p
   
3. s=0,x=0,
   2 W* |( G( Y3 m% d9 H) K, w# ]! X
4. while{x<=1, //while循环算法；
   4 Z! k* c0 C$ e: I9 A, x
5.     y=1, : w$ X% Q- l2 w, E
   
6.     while{y<=2, % F5 l$ |; p9 A7 L# S6 b/ D- ]5 K/ M
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), # L( X6 T+ f# r
   
8.         y=y+0.0009 : @; o! B' h. R
   
9.     }, : ^/ F7 }5 B  Q$ X: I
   
10.     x=x+0.0011 % _+ k+ ^! S1 m
    
11. }, $ f% q$ f6 `5 l- w
    
12. s;7 _. h0 A& `8 a! t1 h0 ~* ?, `
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
7 V0 O; s& w) f. H3 V4 l1 H0.734   //时间，秒
% A; o! |  @9 R- U2 r6 A  @' ?# |
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
0 r/ w. T# K! X' m) l4 B
3 j" V$ X% S9 O－－－－－－－

" f8 _# p. P+ M4 i* X# E+ iForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   2 ?7 M4 d% M7 D& A
2. t=sys::clock();% ?0 I$ S1 e- x# H
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ' u* z; Y\" O5 P7 _1 g: C) F& e
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];  e: M, o5 ?3 z\" i$ p' S
   
5. [sys::clock()-t]/1000;

复制代码

结果：
% T* o! W! B2 I3 X1008606.64947441
5 Q5 ~: f! u+ i5 z$ H2 {0.719   //时间，秒

forcal

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

/ [' q$ P9 T1 B! Dmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  H+ @) r: D& V* b/ W
   
2. tic  y! U( O/ q\" \3 P) v8 |
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);( A8 y2 ^: d; Q& F% B\" F' D! ]
   
4. sum(sum(arrayfun(f,x,y)))
   8 W, t2 D( Q\" I
5. toc
   0 `1 B\" _' w$ Y( S$ d3 Z
6. + @/ S! b. _$ a
   
7. ans =
   8 `! q: j# x  D) L! q
8. % j\" L* \. T5 O9 H- R' N1 T
   
9.   1.0086e+006
   ( a) K! X. g, _& X; c2 t
10. - F8 n% Y2 K& ?6 s
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];1 c+ b# o! r  a. H# T
   
4. mvar:
   
5. # [2 P. i& F. |. l\\" q\\" S7 V
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. * s# R+ A' i5 L$ u8 w
8. t=clock(),
   
9. 7 Q9 e4 t$ P% m3 w
10. oo{
    
11. 8 [# j7 Y, h4 X6 `7 U
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 6 Z# ]/ u\\" x) b8 Q% R
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. & @3 f. I+ m% w: e\\" @+ E
16. };
    
17. ; D$ R# w9 H3 i$ L+ n. l# v$ X1 Y+ N
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">1 c+ b# o! r  a. H# T</font><br /> mvar:<br /> <span style=\\"display:none\\"># [2 P. i& F. |. l\\" q\\" S7 V</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\\">* s# R+ A' i5 L$ u8 w</span> t=clock(),<br /> <span style=\\"display:none\\">7 Q9 e4 t$ P% m3 w</span> oo{<br /> <span style=\\"display:none\\">8 [# j7 Y, h4 X6 `7 U</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">6 Z# ]/ u\\" x) b8 Q% R</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">& @3 f. I+ m% w: e\\" @+ E</span> };<br /> <span style=\\"display:none\\">; D$ R# w9 H3 i$ L+ n. l# v$ X1 Y+ N</span> [clock()-t]/1000;

结果：
) {8 V9 |8 t* ^5 c9 `' L' w% G* ~1008606.64947441
  |2 g: _. }/ a5 g+ b+ p0.735  秒
* i3 c& g  S1 _4 K* H) |) Q" b! @6 \
1 o' k7 Q+ v  x3 K, p+ K可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
" A' _0 j" h$ E; M1 k6 I$ i2 F1 f* Z7 Z
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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