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

forcal

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

1. //用C++代码描述为：
   # c. b7 d0 E& ^; K2 i
2. s=0.0;
   . l6 p\" ]8 _- Q3 D# q
3. for(x=0.0;x<=1.0;x=x+0.0011)
   \" d1 a8 L% s2 C- D( o) |
4. {7 d3 r: X* |' l; T
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   7 J! H' J5 M: [+ y# i+ |0 l
6.    {9 j5 _5 ?* K5 w, r) h
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   7 _8 i# d0 |8 o: S7 _& z
8.    }
   ' t0 q. N8 O& l- W
9. }

复制代码

Matlab代码：

1. tic8 Q\" ]* U0 R5 {* E1 ~
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   - L/ u/ }- P. e, y& c9 n2 U# z
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))6 f- k: q1 o' W  v9 x
   
4. toc
   - Q# f* N( B) b
5. * C( o) m2 `& z4 H1 Q9 l
   
6. s =
   0 R/ F- f3 T+ K. y7 y
7. 4 `; Y8 O\" C7 Q% W' _0 E
   
8.   1.0086e+006
   7 e3 S1 `\" n* @9 i1 Z& R+ Y
9.   I8 K5 O) W$ Y! Y+ V9 N  d
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 1 V5 }8 c% x+ i
4. mvar:6 a% l/ C7 k) }$ @. x; ?/ ?
   
6. t=clock(),\\" \- ^9 [8 l4 G- P
   
8. oo{
   
9. + V\\" T  s8 a2 J* `3 L3 Y* ?5 z1 f
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 N. S  _+ G- d1 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]8 ~; l$ F3 j+ E! P# P6 c
    
14. };
    
15. . F8 A9 l9 i$ E% O' G2 P
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">1 V5 }8 c% x+ i</span> mvar:<font class=\\"jammer\\">6 a% l/ C7 k) }$ @. x; ?/ ?</font><br /> t=clock(),<font class=\\"jammer\\">\\" \- ^9 [8 l4 G- P</font><br /> oo{<br /> <span style=\\"display:none\\">+ V\\" T  s8 a2 J* `3 L3 Y* ?5 z1 f</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">1 N. S  _+ G- d1 Q\\" {\\" ]</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]<font class=\\"jammer\\">8 ~; l$ F3 j+ E! P# P6 c</font><br /> };<br /> <span style=\\"display:none\\">. F8 A9 l9 i$ E% O' G2 P</span> [clock()-t]/1000;

结果：
+ b8 `3 i# p9 l! a" Y1 m+ M" ^1008606.64947441
8 E8 S3 `9 T$ B% G+ ]# D$ w2 ~0.641

* i+ t; Q" L  t5 x) d* M0 oForcal比Matlab稍慢些。
* }/ F1 y/ }; A
－－－－－－－－－－
* ~8 S" O9 [6 P/ {* {
; S. h8 [1 r( Q. K5 y$ Z4 Z再看循环效率。

Matlab代码：

1. tic
   - b6 D/ V9 s- C( I+ A$ w& x7 Q
2. s=0;- T4 W+ \- h5 v0 Z. ?0 u5 d& j
   
3. x = 0;
   9 y! r) T\" {% h) a& t) n
4. y = 1;( z, C& g5 ~. _\" W! q) h& f4 c8 t
   
5. while x<1   
   # D/ C8 v% G' b# V3 |  p) H9 ]
6.     while y<2;        $ b9 }3 C- t; [/ N2 [
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 c' G. i) b5 U8 u) }6 r  J\" V
8.         y = y+0.0009;        
   2 N9 k4 F# V! T1 I& a8 j\" y
9.     end
   8 r: y+ S6 k! ?2 t5 w% U4 w# t
10.     x = x+0.0011;
    8 t: q1 m' u- g$ A6 {
11.     y = 1;
    + M9 ?7 E* a1 e  ^\" o
12. end
    8 k& h3 i( M1 K; k5 O
13. s1 y\" P: T4 z0 C& r9 r' C9 l: L0 J
    
14. toc
    ! _. P' j! u+ H5 b
15. : `  c+ k+ W* W  O
    
16. s =* `4 S: i5 H# b% C2 ^
    
17. ( |) W1 o2 l3 r3 v) p8 u
    
18.   1.0086e+0068 ^8 H4 p, L6 F' A$ P/ O  z. A
    
19. 
    / f% N6 w\" |  T* m' g$ B( }
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   $ F: r, G1 P8 C# q7 y
2. t=sys::clock();* S/ |, w, j( J2 W
   
3. s=0,x=0, 9 P% g8 r7 o) Y9 [9 N0 @0 e- j3 Q$ d
   
4. while{x<=1, //while循环算法；
   & s. ^4 j6 |* g. J+ Q) W
5.     y=1,
   4 l6 t; f6 H2 Z+ X7 t  b
6.     while{y<=2,
   7 S* O; l3 A1 }  n1 J8 B1 g. W' W
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 4 J0 X1 F, `4 k7 d$ c( n6 H
   
8.         y=y+0.0009
   / b& U- Y: ^4 S9 A! X' I8 M5 O
9.     },
   7 a. p6 ^9 O) w( J5 o9 P: [( [
10.     x=x+0.0011
    \" R* A- @1 t+ s6 H* o0 X8 f
11. }, % b0 @! J# G9 C/ D7 n
    
12. s;4 z) g/ `7 R* O0 b, j
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒
. _: X! \4 R0 o  M$ |) U( x
我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
: q: L; n8 R5 m6 u
－－－－－－－

1 W4 y/ G* D# B8 [1 N# T9 ?Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   9 A2 ^( u0 g. X
2. t=sys::clock();, y- D; P% _# O9 N' Y) V
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   3 }% \$ `7 Z9 j
4. sum["f",0,1,0.0011 : 1,2,0.0009];( U. l8 F\" C5 f. q8 V5 U
   
5. [sys::clock()-t]/1000;

复制代码

结果：
) S$ N8 k+ E* r+ A, H+ R( p1008606.64947441
0 {+ I1 x. m: ?3 @0.719   //时间，秒

forcal

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

9 e: ~0 q' D5 i. K/ l; Xmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));: r5 O7 i- _8 ?* s- X6 n
   
2. tic
   9 T. t+ B, S8 W% _. O) }) s, }! l
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);1 Q$ X$ `3 h0 l- c( D
   
4. sum(sum(arrayfun(f,x,y)))4 n: `) m$ J- @& D  i
   
5. toc
   0 _: i- _$ ]$ J) M
6. 
   & ^5 i* \0 q1 d/ I
7. ans =
   4 d( x- O+ g7 M. _! j* y7 `
8. 
   ! b* ~# k' k* f3 m4 V
9.   1.0086e+0068 c/ W( h5 l; i1 y/ H. s( s
   
10. 
    + M( ~; w6 I' t6 u3 Q( @# k  R+ F
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. # [' }! s/ f& ~) H! ?
4. mvar:5 j) s- ~$ B0 B6 M& g3 m3 u6 G
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 0 `7 [$ p1 }1 G6 j  G7 r* _
   
8. t=clock(),
   
9.   B4 U. j, M\\" u7 a, T6 L
10. oo{2 j% h# V3 j0 d& x. k# ]$ p& [( h
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. 1 b, U# I/ B9 _
14.     arrayfun[HFor("f"),x,y].Sum[]# s# i+ W3 s9 Y: o: K& J
    
16. };
    
17. ) @8 d3 ]& f; X: K8 y8 \3 @7 n
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\"># [' }! s/ f& ~) H! ?</span> mvar:<font class=\\"jammer\\">5 j) s- ~$ B0 B6 M& g3 m3 u6 G</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">0 `7 [$ p1 }1 G6 j  G7 r* _</font><br /> t=clock(),<br /> <span style=\\"display:none\\">  B4 U. j, M\\" u7 a, T6 L</span> oo{<font class=\\"jammer\\">2 j% h# V3 j0 d& x. k# ]$ p& [( h</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">1 b, U# I/ B9 _</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\"># s# i+ W3 s9 Y: o: K& J</font><br /> };<br /> <span style=\\"display:none\\">) @8 d3 ]& f; X: K8 y8 \3 @7 n</span> [clock()-t]/1000;

结果：
# O$ Z2 {4 H0 ]1008606.64947441
' U8 X, P8 [* E; d; W+ p8 m% d0.735  秒

1 i6 T5 z6 z/ n, F可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

5 t2 u6 [" X) B: ~－－－－－－－－
# u# S5 K" j: d$ ~$ g0 R+ R2 X# Z
' a" N4 I- d( d) U, j* L1 f, c从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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