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

forcal

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

1. //用C++代码描述为：
   4 \( t; ?9 R1 l$ v3 s' U3 O# I
2. s=0.0;
   6 h+ r* p( j  h9 d+ B, l
3. for(x=0.0;x<=1.0;x=x+0.0011)
   ' V% r  l  N4 }
4. {
   $ P. G$ G, T5 M1 c8 s' @6 [
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   . w2 H$ h8 q/ z/ x8 s6 |% _
6.    {
   7 W  E/ A# G9 k6 L% p0 Y  @
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   : N2 r! k( j3 ?, i: ?* c
8.    }) }3 G/ U  Q* P6 p; [* x
   
9. }

复制代码

Matlab代码：

1. tic9 u; P: M\" i* v8 q% B
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ! S2 f. ?4 N0 ^) o0 ]( _
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   / F& p5 `/ O  `
4. toc
   & H4 ]$ U( n+ C  v. y0 Y\" g- q
5. 
   3 G: i. r/ t3 l
6. s =6 Z) @5 k2 e! S. T# L2 f\" p# x9 t
   
7. ) A5 K( J& h\" Y8 h9 D8 e) d5 @, C
   
8.   1.0086e+006
   . |\" X- _( y, J4 W
9. 
   \" [# Q. N2 ]' e) {
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. + M, w' C, V9 h) V
4. mvar:
   
5. 5 q+ t4 E/ E4 q) r' I) s  N+ u
6. t=clock(),
   
7. . h' M! D2 F+ B$ m
8. oo{
   
9. ) Q: I; i* a; I$ p# z( \' Y
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 R/ j2 l8 t\\" ^1 c
    
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. \\" V6 O* z, t9 N5 W) z0 B
14. };# _3 N2 b; }, O+ b8 s; D! X5 a* v2 C
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">+ M, w' C, V9 h) V</span> mvar:<br /> <span style=\\"display:none\\">5 q+ t4 E/ E4 q) r' I) s  N+ u</span> t=clock(),<br /> <span style=\\"display:none\\">. h' M! D2 F+ B$ m</span> oo{<br /> <span style=\\"display:none\\">) Q: I; i* a; I$ p# z( \' Y</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">1 R/ j2 l8 t\\" ^1 c</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\\">\\" V6 O* z, t9 N5 W) z0 B</span> };<font class=\\"jammer\\"># _3 N2 b; }, O+ b8 s; D! X5 a* v2 C</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
' n# }8 c* [# r2 X; U! Z( ^  v3 ^0.641

* U& F, U9 C) R! q/ j& G7 S) I: XForcal比Matlab稍慢些。

* P7 s, u: @2 c6 @! j3 ]) n7 Q－－－－－－－－－－
6 X1 ]* a1 E" I' H7 a
再看循环效率。
5 ?$ T( h# c1 U  m
Matlab代码：

1. tic
   + y& l! Q) p- M& F' o
2. s=0;; `' N+ T+ {8 w
   
3. x = 0;- {0 I$ ^0 E5 l
   
4. y = 1;
   % K1 m: P' u9 ^
5. while x<1   
   6 Z8 @) ~2 }  c0 I) n2 D
6.     while y<2;        & p\" A% R5 T0 w4 i( Z! o
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));! i1 e% \% F\" e* I% f
   
8.         y = y+0.0009;        
   , c4 ?& X7 [) ?
9.     end
   - }1 E* Q' _% k0 ?, h4 E
10.     x = x+0.0011;
    % L% p2 y) P8 P0 t4 T* P
11.     y = 1;
    1 T/ m\" W' Q$ E, s! e
12. end
    , D  q6 h9 P3 B
13. s7 c0 E4 m4 g0 M; `\" u- c1 U
    
14. toc
    3 }2 ~; M* R6 w- g( D. g
15. 
    / Q\" ^3 j0 U/ }
16. s =
    & g& d! j! ]- F\" @# L; w' [; v0 Q
17. , `+ k% a9 K+ V/ [5 f3 q
    
18.   1.0086e+006
    6 Q/ Z& ~7 d% p( H, T) U9 ^
19. 2 H1 Z% }8 M# y3 e8 t
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   3 ^' W7 j4 g9 e0 g' o7 ]
2. t=sys::clock();
   . O) c3 P+ H0 B$ S6 l. E
3. s=0,x=0, ) V1 X' K4 j9 E, e7 n
   
4. while{x<=1, //while循环算法； ) z: _  N\" e7 F/ U; ]! ~# ]2 b0 K! ^
   
5.     y=1,
   + Y3 Q8 |# ^) S6 X1 d! U; }
6.     while{y<=2,
   / T' @* ^. [; T: Q; `
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   * `' {, A0 k* C9 v. F# Z8 Z( Z
8.         y=y+0.0009
     @; L8 [3 Z4 N) p\" v
9.     },
   3 V5 l7 p% G- O$ _. V% R/ `- `
10.     x=x+0.0011 1 U0 ~4 m& o, e+ I! A$ T2 b
    
11. },
      S# M/ N; i$ B' K) U0 Q; X$ ]\" J
12. s;. D, }2 V1 H- J0 L
    
13. [sys::clock()-t]/1000;

复制代码

结果：
+ x/ e: P: t* a$ q( L. [1008606.64947441
0.734   //时间，秒
- E, |- {! K+ W
/ |2 A* ^- J3 O: m" ]7 l, S7 a8 |我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
- o. p& Z* T! b) W" \+ ?
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:/ M& W\" K$ G7 |% m1 d: S
   
2. t=sys::clock();
   4 g9 E; \5 U4 {
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 5 i* i/ r0 F* \8 F\" K7 ?
   
4. sum["f",0,1,0.0011 : 1,2,0.0009];\" i. g- i6 S5 y( q
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
. H5 M2 w; H- s7 l: A8 Q) `. Z0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
: b% k" p( p, Y) Q+ N( U. j
matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - j# ?) c1 Z- ]$ x
2. tic9 d6 i$ D: E- }# A) y; |# @
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   $ `0 t% V4 ?2 c+ q7 w2 _$ Z3 K
4. sum(sum(arrayfun(f,x,y)))
   $ i- X$ M: L5 n& ?& N
5. toc; J2 u) T# K\" n$ i- p* L
   
6. 
   0 W3 L6 y( v* x
7. ans =\" g! y: T+ }0 ?6 b# }
   
8. 2 @) U; q, g& X# ^- T* y+ i3 f, R
   
9.   1.0086e+006: P* p6 S0 |1 r
   
10. 9 Z9 d& }6 k* A( z) z$ f
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. / ~% X1 B6 G( {7 g  i% C3 K; {$ e
4. mvar:
   
5. * u' J, e. e. _& N0 P: d
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 9 U2 a6 e' I) ~3 R8 N
8. t=clock(),; r; }1 m  n: o% [% r
   
10. oo{- k; s$ k. n5 A; U6 A\\" Y
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],7 h6 f$ f/ Q7 b/ Z7 L5 W3 K
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. % |3 \6 m' |' H+ `0 ^9 n
16. };6 ?7 A! d' v2 c' `
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">/ ~% X1 B6 G( {7 g  i% C3 K; {$ e</span> mvar:<br /> <span style=\\"display:none\\">* u' J, e. e. _& N0 P: d</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 U2 a6 e' I) ~3 R8 N</span> t=clock(),<font class=\\"jammer\\">; r; }1 m  n: o% [% r</font><br /> oo{<font class=\\"jammer\\">- k; s$ k. n5 A; U6 A\\" Y</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">7 h6 f$ f/ Q7 b/ Z7 L5 W3 K</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">% |3 \6 m' |' H+ `0 ^9 n</span> };<font class=\\"jammer\\">6 ?7 A! d' v2 c' `</font><br /> [clock()-t]/1000;

结果：
$ D$ E+ c4 m8 v9 ^1008606.64947441
0.735  秒

0 `9 W9 u, n! f8 n. Z! v/ h可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
) ?. S- {/ [7 Z/ U( B4 t; V# R+ d
$ S' W" ]( }7 x5 |$ r从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。