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

forcal

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

1. //用C++代码描述为：# K8 o, y4 M0 b% l
   
2. s=0.0;
   5 s8 D& g6 w: b6 M! ]
3. for(x=0.0;x<=1.0;x=x+0.0011)
   & k. V. q& y+ c4 E: M# \2 Z
4. {, X1 Y; M$ @; C* g$ i. m
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)! z8 j# ]7 `7 C8 K& g
   
6.    {
   $ c' q9 S, @9 g- a* Z1 w
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));  b9 L1 f* y0 p& X) ]
   
8.    }: ~/ K  V8 r; H/ x2 x\" j
   
9. }

复制代码

Matlab代码：

1. tic) \$ c; ~2 ^9 [' S4 |5 P1 l
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   6 {2 r\" o. b5 b+ z- @! h
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   1 U( t. p+ O6 `; u# I1 y6 E
4. toc
   5 s7 f\" @! v& d1 @) G/ z: _
5. 
   ; |, h3 @  o: ~
6. s =# q2 h8 k, h$ F% `2 P7 @- R' |& B$ P
   
7. ) N! _! A: D- e2 [) V
   
8.   1.0086e+006  _$ [/ {) C, X
   
9. / u; U) k. [5 G8 P& t2 l
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];# G! Z\\" V4 b  Q+ C) D; y
   
4. mvar:2 u6 u: M9 N+ @! }$ @
   
6. t=clock(),8 b) Y7 P: D2 g+ ^, l
   
8. oo{
   
9. 5 v. I0 F8 [9 V7 i2 \  \7 x* @\\" O
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 6 ~9 o- T5 `( J5 d9 b% q
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. 9 q/ z* o4 \+ a% [
14. };' h- |. }$ g5 M% I! F+ g
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\"># G! Z\\" V4 b  Q+ C) D; y</font><br /> mvar:<font class=\\"jammer\\">2 u6 u: M9 N+ @! }$ @</font><br /> t=clock(),<font class=\\"jammer\\">8 b) Y7 P: D2 g+ ^, l</font><br /> oo{<br /> <span style=\\"display:none\\">5 v. I0 F8 [9 V7 i2 \  \7 x* @\\" O</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">6 ~9 o- T5 `( J5 d9 b% q</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\\">9 q/ z* o4 \+ a% [</span> };<font class=\\"jammer\\">' h- |. }$ g5 M% I! F+ g</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
6 u4 I- b5 y: ~. B5 d0.641
* J+ I2 Y# K0 o& Q$ j
; X, {) c: z" g0 ^Forcal比Matlab稍慢些。

! c' w# h8 a2 {2 k－－－－－－－－－－

' E7 F) D8 o" w, G再看循环效率。

( W3 U9 M+ F3 qMatlab代码：

1. tic# j: r* h# U2 s\" T0 F4 S# D
   
2. s=0;
   6 ?: _5 |1 r6 m) k2 b+ g0 `* o
3. x = 0;
   , B7 Q: {/ T- x# L
4. y = 1;# Z1 W1 k9 o/ x4 l( |7 x0 i
   
5. while x<1    & C# R3 i1 ?$ L. Q% K7 F
   
6.     while y<2;        
   6 }2 g$ K& H  E( I8 f
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 S3 |1 x/ s- I/ X. j! Q) F6 \
8.         y = y+0.0009;        
   % j9 i) U7 a% ^: D3 F
9.     end  ^; U; ]\" k; K
   
10.     x = x+0.0011;6 j  U! k: f) P2 b# ~& |& G! A
    
11.     y = 1;
    2 I5 l- P3 C& y2 E) W
12. end\" h. {/ c7 }* h& I) @
    
13. s
    / N6 _$ c  L: R1 |; q
14. toc5 c& H9 _% G, D! p0 z3 ?0 v9 M; r
    
15. 5 O* E* ^, d4 V! o
    
16. s =
    % f) _9 k4 t( ?6 q$ i2 A
17. 
    \" W* R+ o8 h# u, Z4 `
18.   1.0086e+0063 `2 \8 z4 n' @5 {# P1 `: m
    
19. 
    . F0 q6 [\" W9 O( K: o# B
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   $ C$ V, D$ T! T# l+ m) r% h3 K
2. t=sys::clock();# Q  ^' U( S: G- @2 x
   
3. s=0,x=0, ( Q  ^& x# f2 ?: @
   
4. while{x<=1, //while循环算法； 0 R7 B; w\" \4 E
   
5.     y=1,
   % i* X3 p4 c$ e6 S  t* A
6.     while{y<=2, , I/ D4 g2 V( ^/ F; j
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), # l/ |. w0 w! {/ N: q9 m( d
   
8.         y=y+0.0009
   7 Y$ _0 [% p7 T( ?$ w) |
9.     }, ) d. k: U9 Z( S
   
10.     x=x+0.0011
    ; @1 U6 V: Q( D4 f/ {- E) ]+ f
11. },
    ; N2 l1 S( p\" z\" B0 ]. J
12. s;
    ) Y0 l! p! U% ?! _6 u; Y8 E
13. [sys::clock()-t]/1000;

复制代码

结果：
5 v" K$ \2 i( [3 O( x1 ^" V7 n- ^1008606.64947441
8 X! p2 R, ^. ^4 ^6 @' L4 I; Y0.734   //时间，秒

* `) h- j, ^3 {我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
3 h; G$ l6 H9 n. Y+ }3 L
－－－－－－－
/ k8 R7 L" I& [/ Y( q
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   8 i6 Q, q4 x2 C- N' D: S
2. t=sys::clock();' Y& q\" p7 o+ |\" ?
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   / @# p\" o, f: s& Y
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   ( K; V' R; |0 W9 g) \
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
* M' I6 d8 W# t% S/ e6 k3 z+ s; v0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   + t' Y* x6 p3 @4 L! L. E. ]0 @
2. tic
   % Y4 q4 n0 ^# ^% B. ?: n% H5 H
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   1 @, d8 a' y' u/ L3 N$ X
4. sum(sum(arrayfun(f,x,y))); y' \: W) f7 I0 n3 Q5 {
   
5. toc4 }7 ~9 i; i; C9 J
   
6. \" ~, J; T# n/ s2 y\" B
   
7. ans =\" E; t4 R. @+ ^* M! Z- t
   
8. 
   9 J! C; K; C* l* j7 |- y
9.   1.0086e+0069 D* t& Q8 _' s) e\" p8 f! t
   
10. 4 Y8 V2 @$ b% e$ b0 ^0 ]2 q9 o
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 7 d( l% E' h* C\\" M3 z* m
4. mvar:
   
5. 5 M\\" s0 P: v/ O% N+ x4 v/ ~
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 9 a8 n' R5 V, G4 t: S
8. t=clock(),
   
9. 6 y0 x' w* f& q( Y9 B& }8 O9 O
10. oo{7 z/ P+ w/ r! R! N  f8 R
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. ; P\\" o/ c3 I4 N6 b! m) G
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. & [+ }& O- O3 v3 T$ m+ w! v, Z
16. };
    
17. 2 d7 n' `$ S7 x, o
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">7 d( l% E' h* C\\" M3 z* m</span> mvar:<br /> <span style=\\"display:none\\">5 M\\" s0 P: v/ O% N+ x4 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\\">9 a8 n' R5 V, G4 t: S</span> t=clock(),<br /> <span style=\\"display:none\\">6 y0 x' w* f& q( Y9 B& }8 O9 O</span> oo{<font class=\\"jammer\\">7 z/ P+ w/ r! R! N  f8 R</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">; P\\" o/ c3 I4 N6 b! m) G</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">& [+ }& O- O3 v3 T$ m+ w! v, Z</span> };<br /> <span style=\\"display:none\\">2 d7 n' `$ S7 x, o</span> [clock()-t]/1000;

结果：
$ c- Q% }" M. h# l. u; L$ @. V1008606.64947441
0.735  秒

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

. ~. J( B4 `: g) ]4 A－－－－－－－－

0 R- K( A  S5 U# o, u从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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