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

forcal

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

1. //用C++代码描述为：6 g- a8 S) }9 N5 p. {: O3 J2 q4 d
   
2. s=0.0; , K; N8 U4 f, F) y. L
   
3. for(x=0.0;x<=1.0;x=x+0.0011) # h6 v  N! q\" A- I/ T( z
   
4. {2 x% z5 Y  G! K$ P9 L
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)+ J3 ]; T$ M; O6 N
   
6.    {5 Y: ^2 |2 P+ F0 @
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));+ K2 F+ C3 N- s' v5 [) p8 F
   
8.    }
   2 V$ U1 D/ H6 j2 h+ O# T
9. }

复制代码

Matlab代码：

1. tic
   . C. M2 j( a7 l* k1 t. l
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   % L/ @5 b- O3 h8 X7 l
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   4 K7 E( \+ g0 {& ~9 \/ u
4. toc% Y5 R; P5 t( y1 [
   
5. 0 S0 J- G( s9 N) `5 C. o: G8 J
   
6. s =
   ! l( B; N9 L& {* j1 H: ^1 t
7. 
   $ B0 a' u/ O) u9 y; \: K
8.   1.0086e+006
   4 q7 l4 v0 L2 O3 w
9. - R9 Q2 [; O/ c$ r) J1 z5 N+ ~! W
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 9 P0 X. n5 Q- f; P1 e: I8 [
4. mvar:& `1 R+ V  [9 n
   
6. t=clock(),( N, z. Z, ^, P4 Z, E
   
8. oo{, W, l5 i3 q\\" F
   
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],! a7 H4 R1 J\\" ?\\" j0 L) z. s
    
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 j% o7 X, X, s1 K% M9 b
14. };
    
15. & L8 e; Q7 g( @4 r, L- j
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">9 P0 X. n5 Q- f; P1 e: I8 [</span> mvar:<font class=\\"jammer\\">& `1 R+ V  [9 n</font><br /> t=clock(),<font class=\\"jammer\\">( N, z. Z, ^, P4 Z, E</font><br /> oo{<font class=\\"jammer\\">, W, l5 i3 q\\" F</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">! a7 H4 R1 J\\" ?\\" j0 L) z. s</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 j% o7 X, X, s1 K% M9 b</span> };<br /> <span style=\\"display:none\\">& L8 e; Q7 g( @4 r, L- j</span> [clock()-t]/1000;

结果：
1008606.64947441
4 X, Z3 J7 X6 ^+ ^; @- N  o- c0.641

Forcal比Matlab稍慢些。
. f+ B$ E2 g7 v# c) C# D
# V$ H# x. ~% |0 M* A+ E－－－－－－－－－－
# _# A) y; Z8 s/ n# h
! G. C7 t/ r$ }% p再看循环效率。
# K; ?1 m4 G5 j% T$ o
+ R$ P( N+ N6 V0 Z% GMatlab代码：

1. tic
   $ Z& h$ a& w- V# Y5 W
2. s=0;
   ; B: n* p$ _+ M% e# Y! u! r: b
3. x = 0;
   & J1 y\" t) g( ~, E' t5 F
4. y = 1;2 X$ C+ r1 f, W7 Q. T
   
5. while x<1   
   8 x  s6 K4 \/ e
6.     while y<2;        1 p! t+ `, ~( S2 n; I% F
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ; o0 d. O+ |: j: \
8.         y = y+0.0009;        9 m: F0 I, L\" x
   
9.     end
   0 w5 {5 y1 x1 ^* H/ b
10.     x = x+0.0011;
    1 T2 n+ @1 m, E0 g9 b/ H1 M# [: U
11.     y = 1;! U, {! M) Q6 p$ S
    
12. end
    ; W) \& h# W1 I5 y  }9 a
13. s
    6 G2 |2 p. V# l+ X& A1 V
14. toc
    + K4 B! f) w# T
15. 
    4 H4 w/ a3 J6 r6 n( d
16. s =) n3 i: r4 O- a  s) N
    
17. 
    $ I- V6 r) v3 Q) @' b  M
18.   1.0086e+006
    $ I; y- Q' H) z$ i, _1 w- V
19. 7 b7 Z2 n! g; U) Q8 D8 L
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:( H% k2 F/ U3 K( u
   
2. t=sys::clock();
   ( a) z9 u% b2 }2 g5 q
3. s=0,x=0, 3 |4 X\" T\" Y, J
   
4. while{x<=1, //while循环算法； * `9 O' c& p6 [9 _8 T6 i6 @9 D
   
5.     y=1,
   ' b1 T1 U8 A4 I8 \1 D
6.     while{y<=2,
   3 b: U8 Q8 P8 J( |6 B- r9 K4 V- x
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
     h2 Y4 {2 c* H' ]
8.         y=y+0.0009
   + x* W/ m6 k( D6 I
9.     },
   ( w# ^9 M: X# b( S# `
10.     x=x+0.0011 7 m1 X2 t- C' r
    
11. }, ' [2 ^8 i7 M% w' d! n4 Z2 m
    
12. s;* t: Z' X9 P* t+ x: x9 A3 Y
    
13. [sys::clock()-t]/1000;

复制代码

结果：
% R. a; a' S) ]1008606.64947441
0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

－－－－－－－
1 F3 s1 Z# a: j) V7 n8 g$ a& Q* |7 Z
$ c7 {% Y+ E8 g! V6 JForcal中还有一个函数sum专门进行这种计算：

1. mvar:\" j! s$ U6 x9 P8 i, o9 V
   
2. t=sys::clock();# z4 B- O) R# O1 u\" d( O( z
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   , T' j0 i* L- O
4. sum["f",0,1,0.0011 : 1,2,0.0009];/ o  ]) J9 L! y
   
5. [sys::clock()-t]/1000;

复制代码

结果：
% _- S/ J3 H+ p. b! J1008606.64947441
5 h5 L5 @9 }6 F8 l. t' c0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
5 w# P5 I! f; ~
! j2 A( E3 L3 |matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 K7 x; z6 b, G- i% C
2. tic/ K1 V1 b, C8 b. }5 m
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);6 {  e1 X4 [/ A/ P9 l) A* z
   
4. sum(sum(arrayfun(f,x,y))); J* T1 n6 ?3 B0 n  L# m2 q
   
5. toc
   / [* w  }. e4 v$ ?
6. 
   & K& g/ `% u/ q- L' l- r9 `) \* y
7. ans =) x+ ~' a' R  l! Y* G
   
8. 
   7 _) q9 p' B. d, z
9.   1.0086e+006
   ; B* [  @! @3 d\" k
10. 
    & |- z\" y0 G& U7 d
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];\\" J4 z; W0 t5 |* }\\" t
   
4. mvar:$ p! d\\" _- ?- f
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. * @, Z* y9 m9 t8 K
8. t=clock(),
   
9. 0 H! a# }' S\\" w1 i* m& y
10. oo{
    
11. % b% t' n5 Q5 A1 r- A\\" w  U2 ?
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13. \\" n5 l) E* W+ M$ B9 }
14.     arrayfun[HFor("f"),x,y].Sum[]7 u( p' I6 P; h* i3 [
    
16. };; |4 s6 M6 ]\\" k/ B2 I3 I$ k! U  p
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">\\" J4 z; W0 t5 |* }\\" t</font><br /> mvar:<font class=\\"jammer\\">$ p! d\\" _- ?- f</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\\">* @, Z* y9 m9 t8 K</span> t=clock(),<br /> <span style=\\"display:none\\">0 H! a# }' S\\" w1 i* m& y</span> oo{<br /> <span style=\\"display:none\\">% b% t' n5 Q5 A1 r- A\\" w  U2 ?</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">\\" n5 l) E* W+ M$ B9 }</span>     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">7 u( p' I6 P; h* i3 [</font><br /> };<font class=\\"jammer\\">; |4 s6 M6 ]\\" k/ B2 I3 I$ k! U  p</font><br /> [clock()-t]/1000;

结果：
: S  Z' d0 W$ f' r2 F1008606.64947441
( }$ Y, d0 F/ n0.735  秒

& M3 ~% A9 U; b6 o可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

8 g: f& I  M( w3 M! ?- z! b% k% P－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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