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

forcal

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

1. //用C++代码描述为：' s2 N6 Z2 q4 y: Y
   
2. s=0.0; ' `9 f; F: o3 X. a9 W
   
3. for(x=0.0;x<=1.0;x=x+0.0011)
   3 i- P  N# S0 l3 G/ o$ a; j( k& J& o
4. {
   7 w4 p\" {/ k8 M5 Q7 R$ _  l' c
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   ' q. D3 u8 S' `. m
6.    {4 H  c' O, A' L9 K+ V\" [# D& o9 m
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));) L9 r. ^6 N( e& B5 d
   
8.    }
   7 }; u5 Q4 l) o- C9 F# E- ?# ]+ e7 L
9. }

复制代码

Matlab代码：

1. tic
   \" e* y- m7 K% L. p2 d
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   / {  }) N3 q; W
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   + a3 F, [. W' s4 ?& J
4. toc
   : q+ t# i. G1 d4 s8 j) Y2 Q
5. % [( ^1 u; n1 U\" }; Q. \
   
6. s =) V5 O: ?4 r0 ?) a9 U- ?. T
   
7. 
   % f6 \) p6 K( b4 f4 N6 h3 r
8.   1.0086e+006
   8 n/ o4 K2 ^% p% t
9. 
   1 i' `1 |1 p' I) m% t
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];6 J\\" @! G. [% f* d# W/ z; z5 p6 t
   
4. mvar:. y3 J# y: F' r$ |
   
6. t=clock(),/ P+ J( `& m+ Y% R7 O+ R
   
8. oo{+ g7 J+ ~1 b  N9 ?, J: D; R9 Y* O4 T
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],9 d( s1 x2 m1 B, T  V; K( s  ?! p) `
    
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. : h8 P6 q; ^' k: j
14. };
    
15. 3 ~1 k1 h+ N6 T# n
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">6 J\\" @! G. [% f* d# W/ z; z5 p6 t</font><br /> mvar:<font class=\\"jammer\\">. y3 J# y: F' r$ |</font><br /> t=clock(),<font class=\\"jammer\\">/ P+ J( `& m+ Y% R7 O+ R</font><br /> oo{<font class=\\"jammer\\">+ g7 J+ ~1 b  N9 ?, J: D; R9 Y* O4 T</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">9 d( s1 x2 m1 B, T  V; K( s  ?! p) `</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]<br /> <span style=\\"display:none\\">: h8 P6 q; ^' k: j</span> };<br /> <span style=\\"display:none\\">3 ~1 k1 h+ N6 T# n</span> [clock()-t]/1000;

结果：
6 s* z. B# l6 f6 {+ r1008606.64947441
/ L" j+ `; P! u  l( u' A4 v0.641

Forcal比Matlab稍慢些。

－－－－－－－－－－
2 p1 V4 E. t& p
9 P$ e0 c9 H0 U* M再看循环效率。
$ n% k" Z, P/ g4 l
Matlab代码：

1. tic0 S\" }- \- Z1 o; J
   
2. s=0;
   ( }\" q: V( |: U/ l! }, _% {% N' W1 X+ T
3. x = 0;
   ) F& B/ u7 o\" z6 U) e4 @; U0 Z\" g
4. y = 1;  w( d& T) t6 Q* f
   
5. while x<1   
   ; F, Z( H9 q7 O, a+ _, L3 ~
6.     while y<2;        
   . U# t, @$ l' L2 v
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - o* e+ ]4 Z3 u( _, d
8.         y = y+0.0009;        3 h6 G' V/ g( S* Y
   
9.     end
   $ X: o3 c9 v% y0 E& m
10.     x = x+0.0011;' c( x+ e* f2 I! L. _! O
    
11.     y = 1;0 I# O0 O9 C\" h7 y) v( E2 D, _! h
    
12. end
    9 d8 ?7 v; F- I
13. s% M8 i1 `, W% ^, R
    
14. toc
    9 k' {1 t# _! k. L! u% W
15. 
    + ~/ X. h  s5 b) l
16. s =
    1 r6 R+ ?! b) ^- k8 Q2 w: v
17. 8 _8 P) P1 Y* _( G
    
18.   1.0086e+0065 j/ ?  F3 v% b' l! p1 c! [: e7 x
    
19. 
    0 o  X* x7 K0 \) a  c8 U
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:& W2 L0 x\" A+ S& b2 o
   
2. t=sys::clock();# Z/ E% O' x1 b0 m! K7 p' D
   
3. s=0,x=0,
   , x. D3 d5 k3 s9 T! t+ Q
4. while{x<=1, //while循环算法； ; s. t( m% {+ T5 `
   
5.     y=1, 6 `\" v) h# ^- M
   
6.     while{y<=2,
   % `) w: F4 P+ r
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ; G* B$ \; M& U/ L\" n
8.         y=y+0.0009 9 _9 K9 X# z( `: ~) s# L
   
9.     }, ; s3 E% X# ~! O; t
   
10.     x=x+0.0011 . ?3 {' l) q% [4 T
    
11. },
    1 H# Q/ }2 ?0 x9 K7 }- N
12. s;5 U# T/ x# ?0 x( I8 |9 a
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

0 S2 Q; W6 O" {% J, k我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－

Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
   - m' H\" G8 x* A% t0 M' s0 M( X- M
2. t=sys::clock();
   2 u% u( _. ]4 [9 E' G) s
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & B; `1 Q  \% F1 v. q3 Q2 b( r
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   , c% M/ ^6 P' N: x, {* G( x! J
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   . W4 b/ f! g. Q3 A/ q
2. tic
   * l7 v# r' y/ ]! c
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ; F5 l- H\" O& \  l, o
4. sum(sum(arrayfun(f,x,y)))
   2 _% K1 m, u1 @0 I5 Z# ?* ~
5. toc
   - B7 b- P% S6 J& P9 }8 }
6. . s1 Z! Z' }) w. X
   
7. ans =
   6 c3 h0 x+ w+ u; O; b
8. 1 r# G& M3 T6 R+ m4 @+ y5 b
   
9.   1.0086e+006+ k$ e: S  K  p
   
10. , U- f& a2 \5 ]
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];: M6 b% a+ k- R\\" s, a/ W+ k
   
4. mvar:* E2 t3 J3 d* ^3 X6 {1 _$ V) B% r
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. - w' r, |, _- a/ f9 W5 H
8. t=clock(),
   
9. $ \, c\\" [& L1 m/ c7 m2 a
10. oo{
    
11. 0 T! x7 S6 p2 l5 ?5 B% F, f* X9 r# L# m
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. ! R1 A; p; g0 l# ^7 J- ^2 e
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. & }) ]; E$ g; I1 P2 @
16. };. B\\" K2 e9 T8 D4 z! `  x# B* c9 h
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">: M6 b% a+ k- R\\" s, a/ W+ k</font><br /> mvar:<font class=\\"jammer\\">* E2 t3 J3 d* ^3 X6 {1 _$ V) B% r</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\\">- w' r, |, _- a/ f9 W5 H</span> t=clock(),<br /> <span style=\\"display:none\\">$ \, c\\" [& L1 m/ c7 m2 a</span> oo{<br /> <span style=\\"display:none\\">0 T! x7 S6 p2 l5 ?5 B% F, f* X9 r# L# m</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">! R1 A; p; g0 l# ^7 J- ^2 e</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">& }) ]; E$ g; I1 P2 @</span> };<font class=\\"jammer\\">. B\\" K2 e9 T8 D4 z! `  x# B* c9 h</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
! y- V  S& q/ S7 i9 t/ P; v: [0 d0.735  秒

+ ]5 O6 e+ q  d) Q* P5 f1 [7 q& K可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
$ m, g& B8 l8 r) c
9 m7 u* i4 Z. N+ ?/ g4 q$ T1 P& p; ?－－－－－－－－
: e: z( }" \* p$ i
: q7 y# j5 Z1 C8 t" g  w7 E从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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