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

forcal

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

1. //用C++代码描述为：, h% F3 V2 p6 V9 w6 Q& h
   
2. s=0.0;
   3 b% @3 g0 U8 @- i( h4 K* Z
3. for(x=0.0;x<=1.0;x=x+0.0011)
   1 b% M* f# l& ?& D' s- c& z! Z
4. {
   & h0 Z- h* l; A  T7 R' Y
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   , j; |2 e  a# H
6.    {
   : I6 E; Y. h( y, |  f, v6 f
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' m! _* G$ ?9 N7 @
8.    }3 \1 s( }- _' N$ W% \, y, O( }
   
9. }

复制代码

Matlab代码：

1. tic
   4 x0 `; c7 o- ~( {2 t
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);8 H! {5 ^- V3 \6 S5 q: `
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))# t# V7 z7 H% `, Q4 \
   
4. toc0 |+ E\" X) N& m( e
   
5. ; ?8 q5 t; R5 O* }
   
6. s =
   / E+ e5 T3 J9 u3 D
7. ) M2 U& [- U/ Y. l
   
8.   1.0086e+0065 V* O+ I8 f# l( Y0 B0 J
   
9. 
   : B/ ]& a\" I7 H
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 2 E) c1 U2 Q* ^' c( s7 H# a! j
4. mvar:
   
5. + @1 L9 ]4 z5 B( D% O7 o) G
6. t=clock(),
   
7. 8 q4 @1 k. F8 w+ w4 `) {
8. oo{3 r6 l* S0 r, s6 A3 n( b
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 9 {* L* T. y; z\\" D6 |
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. \\" c. ]+ y! s5 N) K8 i3 g# L1 D- ?  A
14. };7 O* p6 \' H) p0 q' t
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">2 E) c1 U2 Q* ^' c( s7 H# a! j</span> mvar:<br /> <span style=\\"display:none\\">+ @1 L9 ]4 z5 B( D% O7 o) G</span> t=clock(),<br /> <span style=\\"display:none\\">8 q4 @1 k. F8 w+ w4 `) {</span> oo{<font class=\\"jammer\\">3 r6 l* S0 r, s6 A3 n( b</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">9 {* L* T. y; z\\" D6 |</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\\">\\" c. ]+ y! s5 N) K8 i3 g# L1 D- ?  A</span> };<font class=\\"jammer\\">7 O* p6 \' H) p0 q' t</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
! w4 _4 w- D& b& _! s7 i0.641
  Y) i' q" F7 m, \9 o
Forcal比Matlab稍慢些。
* n+ V% k0 L" h/ j# H7 n
－－－－－－－－－－
; _) j( Z( Y& w3 O- R5 \
再看循环效率。
' A0 V: V8 k' \9 Q
Matlab代码：

1. tic
   8 J) _7 r8 }4 Y3 _\" ~0 V
2. s=0;
   9 `! a. ?4 V* T/ U5 g! ^) ^2 _
3. x = 0;
   ! I. h5 m' y. x: {( o
4. y = 1;: C. s! ~6 Y& R/ U& j
   
5. while x<1    / M  r* F\" l2 o6 C, O
   
6.     while y<2;        
   8 h* M8 o# q/ C% r8 [$ E
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   8 d) u* h& j* e. ~
8.         y = y+0.0009;        8 A) Z# C4 m+ i& m! _/ d. F
   
9.     end& H( i+ z% r& t* a: w
   
10.     x = x+0.0011;! Y3 g, f. w+ R& k
    
11.     y = 1;  G# C. [+ L# W2 w5 l) K
    
12. end
    4 }& a% x1 L, h  u\" H+ c
13. s
    1 f: l7 a! Q) U( u9 Y, {% ?7 g
14. toc
    0 r, F- [2 I. z
15. \" ?9 u( p/ o6 Z0 R* L$ n
    
16. s =
    3 C* P, g' o$ M8 n, m. V% e
17. ! T1 G# S! _, }4 W
    
18.   1.0086e+0067 K- o# o0 n/ d5 x' j  q
    
19. 7 f. c: G' a9 C
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:5 ?$ s; V0 O5 Z  x- R/ |9 m2 N
   
2. t=sys::clock();
   * T\" ?- O1 M* p# E% x: w3 n\" E
3. s=0,x=0, 7 Z; y\" E% {% X, `; B
   
4. while{x<=1, //while循环算法；
   \" V: N% i$ Q, `# L0 c
5.     y=1, 2 Q9 A! C/ s/ ^& ]$ L7 z& @8 u9 |* I
   
6.     while{y<=2,
   : q# r. S- V% s/ A' B
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   1 U7 I) c' L  X- T# @
8.         y=y+0.0009 - N) v$ |1 H+ [) l  ^8 o
   
9.     }, $ r: c9 H+ v+ w4 z4 S\" Z
   
10.     x=x+0.0011 1 x9 {3 a7 c\" v! A1 g
    
11. }, \" Y: k! x2 f. T  T. N4 a; R
    
12. s;  ]9 A' ^- \9 T( e! ~! E7 Y! A: Y
    
13. [sys::clock()-t]/1000;

复制代码

结果：
" b) u1 i4 R+ u; C* E1008606.64947441
1 r7 y9 \0 S. j1 @4 D2 {" N0.734   //时间，秒
" P+ }7 I& k- c8 {; @: ]) Y
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
4 B0 H! e4 {% W* K' O/ n- S
+ `4 X/ h$ ^2 f－－－－－－－

1 g/ D& X; E- f! K- C7 A  CForcal中还有一个函数sum专门进行这种计算：

1. mvar:7 w8 W8 w$ j1 N% Y: y
   
2. t=sys::clock();! J( E: x7 Q2 e/ V
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   : Z' S+ P; b2 ?* I7 |5 x
4. sum["f",0,1,0.0011 : 1,2,0.0009];2 g6 F/ h( ^\" n: Z
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

# R+ c* p" g: }0 ]  Jmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));& o+ k0 i; D& y4 i: G) M
   
2. tic! g9 m- Y; A' d4 c7 c
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ( O9 `\" g! G( h0 A3 c\" b
4. sum(sum(arrayfun(f,x,y)))
   . C, h! I( d1 r, f$ _$ ~% i' G
5. toc
   # ?! q- o1 v* q1 n! r
6. 6 c& a9 {* e5 I\" v1 S2 P
   
7. ans =
   1 M! C6 S/ Q% V* Q+ ?
8. , I# T! U6 x! P, B$ s, s
   
9.   1.0086e+006
   \" L# L4 G# I& Y/ Z7 Q$ w
10. 
    ; \$ y8 p5 x7 x1 `2 m
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 5 C2 j( q6 e) x: U6 d8 r2 P; j
4. mvar:\\" d7 [# R+ c, Q+ J: z
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 3 f- n$ ?  w. f, s! F
8. t=clock(),. s6 |' @$ v\\" W# K6 T  X# Q( a\\" f4 e
   
10. oo{
    
11. : B\\" N0 q3 K9 p5 n\\" r! N$ W
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],: I  s' b4 C4 x6 C$ v4 a\\" [
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. $ x1 e0 c0 ^4 n$ B9 L
16. };
    
17. ; w$ Y\\" a  W, K8 E) I( u& s' o
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">5 C2 j( q6 e) x: U6 d8 r2 P; j</span> mvar:<font class=\\"jammer\\">\\" d7 [# R+ c, Q+ J: z</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <br /> <span style=\\"display:none\\">3 f- n$ ?  w. f, s! F</span> t=clock(),<font class=\\"jammer\\">. s6 |' @$ v\\" W# K6 T  X# Q( a\\" f4 e</font><br /> oo{<br /> <span style=\\"display:none\\">: B\\" N0 q3 K9 p5 n\\" r! N$ W</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">: I  s' b4 C4 x6 C$ v4 a\\" [</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">$ x1 e0 c0 ^4 n$ B9 L</span> };<br /> <span style=\\"display:none\\">; w$ Y\\" a  W, K8 E) I( u& s' o</span> [clock()-t]/1000;

结果：
& g! x5 }- P* Y1 i5 F8 A1008606.64947441
# W! L0 C7 Z1 |. Z" r0.735  秒
3 D7 w4 f+ z7 |" P" z7 {
( p) g4 a2 p0 H4 c可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－

, Q: [; T( M! @# `! g从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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