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

forcal

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

1. //用C++代码描述为：
   ! N! {5 R* C( c/ k! O! p& x
2. s=0.0; 1 ]9 J: M+ N* H; x& @
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   7 |7 `* W2 l% A+ B4 T7 j
4. {
   8 O% p: i0 V# Y- d1 J
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   \" L, m, ]$ N& ?4 Z8 W( }8 q
6.    {! `: i: H% m1 W, Q9 N% v' ?% c
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 x; C) b! k8 a8 `- I6 X
8.    }$ T& ]! H) ?6 C4 z
   
9. }

复制代码

Matlab代码：

1. tic
     h2 m  t3 X  M# w
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   5 {7 A5 Z$ K$ n: ?% s
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   \" v( b+ W% \. T% `) N  w
4. toc. K# ]; a3 Z$ ]2 C
   
5. * F! @) j. b- v) Y# N
   
6. s =
   1 ^4 o. m. c0 q; r5 Q3 t
7. : g, g* N2 e( S. ?' {) q
   
8.   1.0086e+006
   - S1 V8 x* G1 N/ V# g5 G. }2 G
9. 1 C& q  P. ]. [2 l
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ; i+ a0 h+ I. K% C% B
4. mvar:& r/ Z8 K4 V+ X5 s6 z; `8 x
   
6. t=clock(),
   
7. 7 O7 l& i- Q5 w3 n
8. oo{& D\\" ]# _: J0 ^4 }: j  |$ d
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 1 g4 D' \: E9 H  d# `5 p6 N
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. # ?( \. q& w1 W8 I\\" p; a
14. };
    
15. : `3 p7 B0 |3 G! F
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">; i+ a0 h+ I. K% C% B</span> mvar:<font class=\\"jammer\\">& r/ Z8 K4 V+ X5 s6 z; `8 x</font><br /> t=clock(),<br /> <span style=\\"display:none\\">7 O7 l& i- Q5 w3 n</span> oo{<font class=\\"jammer\\">& D\\" ]# _: J0 ^4 }: j  |$ d</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">1 g4 D' \: E9 H  d# `5 p6 N</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\\"># ?( \. q& w1 W8 I\\" p; a</span> };<br /> <span style=\\"display:none\\">: `3 p7 B0 |3 G! F</span> [clock()-t]/1000;

结果：
1008606.64947441
0.641
! ^* J8 S* c0 m9 ^% b  n- c
Forcal比Matlab稍慢些。

' ~- D, j- C8 I5 `/ j/ B－－－－－－－－－－

再看循环效率。

/ K+ ~, y# \  q6 ^Matlab代码：

1. tic; x3 l' I. P9 s
   
2. s=0;7 E  V3 k) I# X- f2 J
   
3. x = 0;
   9 A  w4 [) c1 @# e& Z* ~) Z; \
4. y = 1;2 q\" ^; k0 @8 E1 d\" f9 I
   
5. while x<1    . H' }- E1 s+ ?4 C& }, p
   
6.     while y<2;        1 H6 a' K& J' x0 z. d( l7 U0 {
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));2 e! y+ i\" z7 s9 g0 b7 n$ a! p
   
8.         y = y+0.0009;        . p  L5 \- S3 h% Q# v/ X! i% r- z
   
9.     end
   4 K: p\" k  G( s% I! }
10.     x = x+0.0011;
    8 h5 A: H6 G/ |  J+ x' U  J
11.     y = 1;! h! Y' {% I% H
    
12. end7 s+ H\" W% H7 W) `. @! J- j4 H
    
13. s
      E9 M\" ~7 Y! i' S- ~6 L
14. toc  T/ i( g0 o6 H, b9 `4 L
    
15. - t6 k9 m6 D: Q2 v4 U! d
    
16. s =
    : v4 ~4 V4 H7 X# m5 K( u2 b9 _
17. 
    ; @# c1 `( D/ f& K* A& u& |; k
18.   1.0086e+006
    9 n! \7 a1 C/ A+ N) |
19. - X\" B( d4 b1 Z
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:2 e) ~1 s\" v' ~/ r$ \
   
2. t=sys::clock();
   4 T& B& T( L1 l! k
3. s=0,x=0, : c& B* v# M  ~6 [
   
4. while{x<=1, //while循环算法；
   2 @5 g4 c) H* f, w! Y
5.     y=1,
   3 }; N# I: {5 p& m0 U% D
6.     while{y<=2,
   7 l9 Z( K\" v2 C$ f! F; I, i+ e
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   3 j! ^6 N$ u; v/ Y
8.         y=y+0.0009
   $ I3 q5 P# G9 s% [
9.     },
   - I& n\" A1 V% I9 a& J' F
10.     x=x+0.0011
    5 g5 S3 U4 G, ~1 x1 F
11. },
    6 Y: m+ k1 A$ o2 ~* S5 G  t1 Q: T+ L
12. s;
    1 q' x$ D1 H$ P
13. [sys::clock()-t]/1000;

复制代码

结果：
; R; j8 X; A4 r1008606.64947441
8 I% ]3 t- T/ S3 I* g3 y' }! E0.734   //时间，秒
: [8 G1 h5 ?. n' O. W2 l& J* m
' p# a7 C' Y7 M, X* v8 j我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

4 i: B4 \2 \- `/ f－－－－－－－

4 ~1 b* N0 S- y0 m9 t. a0 jForcal中还有一个函数sum专门进行这种计算：

1. mvar:
   # _% @9 H0 u! l+ A( Q( n
2. t=sys::clock();5 V; c. m\" J\" X# h! R- p' e
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 O  K& i2 T+ J$ @5 s6 K
4. sum["f",0,1,0.0011 : 1,2,0.0009];\" C; r% d7 F. r8 Z/ H
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
% G3 U3 C$ g& p0.719   //时间，秒

forcal

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

5 y* Z( b: }! ]! T2 t- m% ^$ o; smatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   0 p* W0 ?- F! j. u' J/ B
2. tic: Y  n. ?5 e  D/ m0 T
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   + [! A6 K( X! u, y8 E+ Z8 }5 H
4. sum(sum(arrayfun(f,x,y)))
   , {5 K* f( Y6 e! t5 _
5. toc; E8 Z1 L2 V! Z1 w
   
6. 
   : R* n( p  C2 h& q
7. ans =
   . Z0 Q. Q, W' l! l9 z3 h% D\" I
8. 
   3 @$ R' w% c  H1 E+ u
9.   1.0086e+006
   5 W. l' k) f7 \- c  ]
10. - X: a' G: i1 A/ c6 z1 M\" q
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];# E  |$ ]# p: F4 w
   
4. mvar:2 A/ ?% Z' O9 m: U1 N
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 4 k* s; i! n$ h4 O5 }
   
8. t=clock(),
   
9. # P9 ~! ?) H, t6 C
10. oo{
    
11. ' R! q% h4 q7 N% N0 S
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],8 L7 V: P3 ~6 G# x* }! b5 K3 b1 G8 ^
    
14.     arrayfun[HFor("f"),x,y].Sum[]& B- u7 g- @! H3 l2 v
    
16. };
    
17. 8 L$ m0 \5 \- A+ A1 y9 _
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\"># E  |$ ]# p: F4 w</font><br /> mvar:<font class=\\"jammer\\">2 A/ ?% Z' O9 m: U1 N</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">4 k* s; i! n$ h4 O5 }</font><br /> t=clock(),<br /> <span style=\\"display:none\\"># P9 ~! ?) H, t6 C</span> oo{<br /> <span style=\\"display:none\\">' R! q% h4 q7 N% N0 S</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">8 L7 V: P3 ~6 G# x* }! b5 K3 b1 G8 ^</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">& B- u7 g- @! H3 l2 v</font><br /> };<br /> <span style=\\"display:none\\">8 L$ m0 \5 \- A+ A1 y9 _</span> [clock()-t]/1000;

结果：
* h0 s; j, t/ |1008606.64947441
6 a' p. \9 B* g) e0.735  秒

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

2 x9 N" J5 X, \& p－－－－－－－－
4 A0 A5 l$ D6 ]/ N) o
从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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