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

forcal

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

1. //用C++代码描述为：
     r- s- P3 u8 U9 y6 s8 P5 q
2. s=0.0;
   5 f! ~( `' Y\" G% [; v
3. for(x=0.0;x<=1.0;x=x+0.0011) ! s2 t, }$ W% r3 Z. u
   
4. {. ^/ `2 n# {( c
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   . D5 b. M! _2 d# w4 ~! |( f
6.    {0 |0 w4 H# A) ~0 |+ {$ E: K
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   6 e& K, M- \1 `) A4 j% {
8.    }& ^7 `\" ~8 ?$ X
   
9. }

复制代码

Matlab代码：

1. tic, n+ O4 W$ I/ H, k
   
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   , y3 G) {; m! ~& m
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))% V7 L& E4 _+ b  T- A% @
   
4. toc& x4 {! F4 ?: O7 V+ m  P  q: G
   
5. 
   8 T# x: `3 d+ w4 F, w
6. s =7 o; l\" @, `, s, [/ {5 S
   
7. 
   9 Q& P+ u6 }- u
8.   1.0086e+006) u! X3 n% U: U9 x( a
   
9. : z0 s. L\" V! r! F8 X
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];5 D# j4 k: D8 Q6 s
   
4. mvar:
   
5. ' L, q7 l8 e6 \$ A
6. t=clock(),2 E' X3 r9 N3 X7 C1 @7 V+ O' I0 L
   
8. oo{
   
9. / L& M/ \2 l* q0 N% n& h% ^2 K5 V  N
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],) T6 U5 C  {# |# K$ M
    
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]% L# `: [3 r& \' ^1 y+ E
    
14. };  `) x2 Y7 T3 I, J& f
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">5 D# j4 k: D8 Q6 s</font><br /> mvar:<br /> <span style=\\"display:none\\">' L, q7 l8 e6 \$ A</span> t=clock(),<font class=\\"jammer\\">2 E' X3 r9 N3 X7 C1 @7 V+ O' I0 L</font><br /> oo{<br /> <span style=\\"display:none\\">/ L& M/ \2 l* q0 N% n& h% ^2 K5 V  N</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">) T6 U5 C  {# |# K$ M</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\\">% L# `: [3 r& \' ^1 y+ E</font><br /> };<font class=\\"jammer\\">  `) x2 Y7 T3 I, J& f</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
8 q, M/ y& b! B2 o0.641

5 A7 L& d  y6 g$ YForcal比Matlab稍慢些。
" m. ?; \; h. v+ k+ U
－－－－－－－－－－
" V7 |# o. w4 T
5 {7 [. s* N2 D( Z* x再看循环效率。
; f: j: p1 S* s5 n
7 b. b1 S: q' K6 |1 e9 t# F+ WMatlab代码：

1. tic
   9 D% Y% I9 }) i* k& t+ H* O. b5 M
2. s=0;
   ) U4 ?0 o- Q# `) V
3. x = 0;% p) n6 z) E! X8 u* t3 R0 V! o
   
4. y = 1;! U, e# \5 \' c  I2 p# q
   
5. while x<1    * v. D. H- w\" I
   
6.     while y<2;        3 l8 M; R4 Z- Z9 i3 a* d
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));7 I9 X' {& o0 y; L8 ~
   
8.         y = y+0.0009;        
   6 R6 ~5 U5 F; k0 I$ ^. q
9.     end
   9 c# u5 r# k\" p7 \5 @; u. ]$ ~
10.     x = x+0.0011;% ~( u' K4 @  i
    
11.     y = 1;
    ' W4 B% s\" E8 \/ D8 C4 @
12. end+ Q+ o# }% }% J: y, L
    
13. s1 }5 ]0 n0 _4 t+ L
    
14. toc
      l, `/ k0 F( j0 g
15. 
    9 k, q) u% C9 U
16. s =
    . G4 r  j6 Q/ _& s4 J
17. 
    7 U% N0 E5 l) l) J* w0 ]
18.   1.0086e+0065 I' @7 Y1 C: }
    
19. ) m! R0 z3 g) k2 K0 [* q9 O& U
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:  |; j( b+ s; ~$ f7 E
   
2. t=sys::clock();7 E& D  _  Y6 t\" Z4 S4 \4 c: ?8 y$ v
   
3. s=0,x=0, ' j/ P* Q$ @- B5 g9 }. Z
   
4. while{x<=1, //while循环算法；
   ! G* z1 `$ p, }1 f
5.     y=1, 6 U: Q9 s4 [0 g$ |$ W. `+ `5 ~0 B
   
6.     while{y<=2,
   7 J' V& f0 q% Z\" M- \4 }( g. f
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), + n! J) M9 l; ~1 ^
   
8.         y=y+0.0009
   1 k4 n, t( T3 _- R9 {1 U
9.     }, # x0 M% A* i/ J! W! J( u# I2 p
   
10.     x=x+0.0011
    ) x& ]0 t- g. e* o  M\" b3 V
11. },
      }2 X, J3 t+ {& r3 V
12. s;, C6 X! u! e6 u7 P1 C) t: ?1 b6 i) ]
    
13. [sys::clock()-t]/1000;

复制代码

结果：
/ c9 }6 \$ k) f; M1008606.64947441
0.734   //时间，秒

" V5 m) N; Y; k  k% W我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？

5 f; E* X8 U( S4 R1 k－－－－－－－
! V+ {4 ?* O$ w" {9 _" n; L
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:/ k/ M( Y! Q: O' ]4 T
   
2. t=sys::clock();
   \" x8 L0 E$ _% _$ B3 d) ]7 N
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   5 [5 D. L$ C9 d( E* F
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   4 A  N- {/ N4 S0 X$ S  V
5. [sys::clock()-t]/1000;

复制代码

结果：
. `5 O  r! d+ R  h- y: O6 c1008606.64947441
( M3 d, D$ N( \0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
+ L: N" h1 O/ q! S0 g) l" c
3 q- O; J; o' t+ dmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 D8 f6 j; ^* w. o9 x, d) S9 z% ~
2. tic
   * i# Q% d$ K, P* x& c6 S
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   3 Q- P1 ?7 L, L* n+ z9 ?
4. sum(sum(arrayfun(f,x,y)))3 y# d  \0 t! N& c- M' _
   
5. toc$ d2 Q5 D* m# A
   
6. & N0 R5 N! }( X
   
7. ans =* q/ O& D' m, l: ]
   
8. \" Q0 V& v. O7 j$ h& o
   
9.   1.0086e+006  Y  |* ]7 F! N0 L; P5 O
   
10. 7 E2 K4 S9 Q. n2 d
    
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. ( k/ C\\" i( H; G4 ?
4. mvar:
   
5. \\" g1 C, r5 o  D2 d9 B
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   
7. 1 x* e. n3 ^2 K% ~
8. t=clock(),0 s\\" m\\" j# z7 W5 D
   
10. oo{
    
11. 1 _2 Y! ?# g1 ^$ b. ?! A
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],9 O# M; D; }/ T- q2 c9 p- _
    
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. + v7 w- {+ n- d
16. };* ?3 |, a, Q7 U( t
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">( k/ C\\" i( H; G4 ?</span> mvar:<br /> <span style=\\"display:none\\">\\" g1 C, r5 o  D2 d9 B</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\\">1 x* e. n3 ^2 K% ~</span> t=clock(),<font class=\\"jammer\\">0 s\\" m\\" j# z7 W5 D</font><br /> oo{<br /> <span style=\\"display:none\\">1 _2 Y! ?# g1 ^$ b. ?! A</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">9 O# M; D; }/ T- q2 c9 p- _</font><br />     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\">+ v7 w- {+ n- d</span> };<font class=\\"jammer\\">* ?3 |, a, Q7 U( t</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
. a- r5 x2 ?$ M- T! f0.735  秒

) }% h4 e/ v' ~  X5 O: ?% M( M可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
0 w, x* ^$ O! Q1 H4 f' l; j! h! d
－－－－－－－－
$ V! \- R7 z# x/ n2 ?/ P" W% @
& w3 p# x9 o2 m从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。