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

forcal

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

1. //用C++代码描述为：3 B; z+ x( g* Y9 m7 ~\" g
   
2. s=0.0; 7 L+ o! `* O  @7 g' X3 ^& A# ?
   
3. for(x=0.0;x<=1.0;x=x+0.0011) ' a1 ]+ d- _! z% ?; ?* ~. W
   
4. {
   ( S7 k9 [8 i& A; V1 z
5.    for(y=1.0;y<=2.0;y=y+0.0009)  Y2 \& _( r/ ?8 t$ _5 s
   
6.    {
   \" _# g5 t, E# h9 Y
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   - ^5 a0 `( J6 t4 i& X
8.    }
   , E: J+ J% v% ^6 @, j% \
9. }

复制代码

Matlab代码：

1. tic
   ; t) U- f/ V. ^8 _& y; }1 I
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, X7 S) v+ h$ ~+ h: w6 b- J/ w# W
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   . N% U- G0 Y& X. b/ u! F
4. toc* F8 B8 S8 M\" D7 a1 b
   
5. 2 J3 O7 y: d+ F! Y& @
   
6. s =6 X$ ^' H1 V: O' h  k. m4 V( `& t7 ]4 N
   
7. 
   9 \. M( Q7 J# B7 g1 Z4 E
8.   1.0086e+0067 U0 k& A6 n6 u7 ~* g6 ]
   
9. 
   / I* j\" g0 p# Y8 l
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 6 z% n/ ], F6 Z\\" ~
4. mvar:
   
5. * Z. a5 d  U: [7 V$ {
6. t=clock(),5 }7 w4 [4 N: F0 m  d
   
8. oo{: Q% `( d& H$ K3 H/ d
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],4 D$ z, }( l& ^+ t/ n\\" O3 k+ H. ~
    
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 g7 T/ A\\" Y2 F$ k
14. };  S( H1 ?$ a) F( d3 N
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">6 z% n/ ], F6 Z\\" ~</span> mvar:<br /> <span style=\\"display:none\\">* Z. a5 d  U: [7 V$ {</span> t=clock(),<font class=\\"jammer\\">5 }7 w4 [4 N: F0 m  d</font><br /> oo{<font class=\\"jammer\\">: Q% `( d& H$ K3 H/ d</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">4 D$ z, }( l& ^+ t/ n\\" O3 k+ H. ~</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\\">. }9 g7 T/ A\\" Y2 F$ k</span> };<font class=\\"jammer\\">  S( H1 ?$ a) F( d3 N</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.641
. v( r# Q6 v1 G% a% y& _1 `, b
/ D+ u. e/ d2 ~$ I* N& T* Y; jForcal比Matlab稍慢些。
& c# U: q1 k5 t  A6 G3 T
& h& P9 K- G5 F' r# C2 Z+ D－－－－－－－－－－

再看循环效率。
+ ]- [& z) a- c  h8 A& \
6 c% E- @* b! PMatlab代码：

1. tic* O. p' I' B2 |+ p1 g% @
   
2. s=0;
   8 F( q) I3 H$ l\" G( K% {  K
3. x = 0;
   - G% k: g/ `. Y  d
4. y = 1;
   + W  ~1 q* p+ b5 o% p9 |' h
5. while x<1   
   * @$ v) R# e% r+ `0 k4 ]
6.     while y<2;        ! C6 @* Y4 e5 ?6 e$ t6 r
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' Q5 Y) t+ ~6 q1 H
8.         y = y+0.0009;        
   : X\" g# N* W2 B$ |- U1 a
9.     end
   . v- ^; B) w3 h
10.     x = x+0.0011;
    8 [8 P+ Z2 z9 ^: q
11.     y = 1;
    ; H  C8 h: Q4 ]8 i3 E: z5 S
12. end
    ; e9 U% M% i\" M5 w
13. s7 |! j\" `\" j; a, t
    
14. toc
      K  y6 n. ]( {$ D9 @( x
15. 4 r( R8 o& N6 H
    
16. s =8 C\" u! C. A% ^5 L/ ~
    
17. 
    / C1 X4 p9 d4 w; _/ y
18.   1.0086e+006
    ( ?8 U4 \# E3 {$ P5 ?  _: O
19. 
    ' ]$ U7 |' ~* Q; u
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:& X; [  @  r7 n0 I9 I
   
2. t=sys::clock();8 R. H! L7 |2 y& A+ H
   
3. s=0,x=0,
   1 k/ ]3 d5 q/ V: S6 u6 S$ R
4. while{x<=1, //while循环算法；
   & `1 z0 x\" D4 R$ b
5.     y=1,
   / P* A+ [+ ^5 A: _  n! @
6.     while{y<=2,
   , P7 L* p8 X- J7 t1 d
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   / l8 j4 |3 M% ~4 C, ?' v
8.         y=y+0.0009
   # [& ^( A5 K6 k$ r
9.     },
   : m6 P9 m+ z1 ]- ?2 b7 [
10.     x=x+0.0011
    9 A/ B2 L- h! O) L7 Q' A% v0 \
11. },
    / z0 ]6 j' J' C, k: V( ]$ q
12. s;
    ! K# M) J$ u, i0 q, A
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734   //时间，秒

% K; \3 w+ V8 A0 @1 `我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
; r2 J+ o3 `9 z1 S# r( I! r/ U( O
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:! H9 k% \( E0 v! y\" y
   
2. t=sys::clock();: f+ N+ P9 z7 Z& R0 G& n! ]
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   \" A, ^! W. a7 K1 I% ~' X+ p
4. sum["f",0,1,0.0011 : 1,2,0.0009];+ s6 {. x$ o  `& G8 @5 @& X
   
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
6 K% h( A+ g9 |7 ^
  C$ l9 H& d5 q, i5 nmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' t: i; ~\" L; f1 X
2. tic2 W$ P& U- e0 g% g& d1 L0 g
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);. b+ }8 B' ]2 M  c
   
4. sum(sum(arrayfun(f,x,y)))5 O2 c  G' o\" ~! l* A\" Z; y
   
5. toc- K! a\" A& t# T( q& w7 v
   
6. ( ]/ }) Y9 y4 {8 B! p' \4 `
   
7. ans =/ O/ ?, m\" b6 {( ]
   
8. ( X9 [2 B: E' g
   
9.   1.0086e+006
   # L2 _6 V) e# s  j2 M6 N: I
10. \" x7 T- g2 M! X! q5 ~5 N9 p( i8 J% x
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. \\" o! |  O- u1 o; x; G6 f
4. mvar:+ ~. l) Y7 S4 w) ]' ^
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); * S* i4 ~0 B  i
   
8. t=clock(),
   
9. # N  s% g6 {7 G  u0 o, o! a
10. oo{
    
11. ! [) A( U, V% `; g! a6 {6 T' E6 P- o9 K' p
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. . ^, ]) B# N8 S2 u, `/ m
14.     arrayfun[HFor("f"),x,y].Sum[]; d9 `/ P4 V5 `( N% m& ~6 z
    
16. };
    
17. , _8 t* p8 z+ I1 D, @3 o! c
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">\\" o! |  O- u1 o; x; G6 f</span> mvar:<font class=\\"jammer\\">+ ~. l) Y7 S4 w) ]' ^</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">* S* i4 ~0 B  i</font><br /> t=clock(),<br /> <span style=\\"display:none\\"># N  s% g6 {7 G  u0 o, o! a</span> oo{<br /> <span style=\\"display:none\\">! [) A( U, V% `; g! a6 {6 T' E6 P- o9 K' p</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">. ^, ]) B# N8 S2 u, `/ m</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">; d9 `/ P4 V5 `( N% m& ~6 z</font><br /> };<br /> <span style=\\"display:none\\">, _8 t* p8 z+ I1 D, @3 o! c</span> [clock()-t]/1000;

结果：
$ Q1 z& g. \9 t1008606.64947441
0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
2 b* |& t0 |0 @! u/ K
6 J6 w: ?& Y! u$ L( j: E; K- y－－－－－－－－

( j  T, M' n3 r从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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