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

forcal

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

1. //用C++代码描述为：
   / u- B1 X7 V# J2 ^\" ?
2. s=0.0;
   5 g- P1 }: x! l) d$ ^+ q
3. for(x=0.0;x<=1.0;x=x+0.0011)
   \" t4 n. j4 g' u\" B# t' g3 n8 U* b
4. {
   * `\" o+ j% v9 d5 [6 n8 \2 T
5.    for(y=1.0;y<=2.0;y=y+0.0009)
   - K* b# ^2 T9 u
6.    {
   # L& q! n+ B( Y. B
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & @  u  I$ @\" Y! z\" x; a
8.    }
   ) G\" W5 S* X3 l1 _# P8 L- f
9. }

复制代码

Matlab代码：

1. tic
   + T) r& c2 p0 T( U1 ]. R  P0 ~
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);* x; B. B5 B2 z0 N
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   7 L1 ]  T3 ?3 M: s
4. toc+ Y  s. q/ K6 P% w. u: `
   
5. $ e8 w9 R- s; m1 J+ X. c  y8 ^
   
6. s =( R9 z6 f4 r2 j( w+ w2 P
   
7. 
   3 }# r- Y  g) G, v+ t9 O
8.   1.0086e+0066 m) g; G  s% h. d- ]2 }- \
   
9. 9 ~( L8 G  _8 X1 f0 g8 j2 L
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. - T( k6 y+ c1 q$ @1 J* K
4. mvar:
   
5. 8 l  u$ I+ d- n  Q+ ^. D; x
6. t=clock(),3 D4 B\\" d: ^\\" W4 P
   
8. oo{
   
9. 1 q# @$ q0 k: w* B! i% \0 j
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 J2 l; H' 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]
    
13. . k9 p7 H3 a' g( e3 u
14. };
    
15. , x! ~2 ]+ Q$ A8 `8 E7 E6 h& O
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">- T( k6 y+ c1 q$ @1 J* K</span> mvar:<br /> <span style=\\"display:none\\">8 l  u$ I+ d- n  Q+ ^. D; x</span> t=clock(),<font class=\\"jammer\\">3 D4 B\\" d: ^\\" W4 P</font><br /> oo{<br /> <span style=\\"display:none\\">1 q# @$ q0 k: w* B! i% \0 j</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">1 J2 l; H' 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]<br /> <span style=\\"display:none\\">. k9 p7 H3 a' g( e3 u</span> };<br /> <span style=\\"display:none\\">, x! ~2 ]+ Q$ A8 `8 E7 E6 h& O</span> [clock()-t]/1000;

结果：
, V" R0 c- \# c! p) i: V1008606.64947441
4 A  @; l! x; C) G1 x0.641

Forcal比Matlab稍慢些。
1 R% ~* x5 o7 c, G; M6 D
; W- g; V; n& `' t－－－－－－－－－－

+ `' c, s" |9 c- Q* f+ M再看循环效率。

Matlab代码：

1. tic% M- A4 W) l, S2 J6 ~8 S
   
2. s=0;, }5 D' z\" q; n
   
3. x = 0;
   7 d7 B# U2 u8 j
4. y = 1;2 j* r# z& n: ]6 `, |
   
5. while x<1   
   3 P0 Y2 [7 I( T$ M3 B& r
6.     while y<2;        
   4 d2 D( d6 _' n
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   # s1 S9 x  Z9 K! P$ n$ ?0 x$ ~2 ^6 A
8.         y = y+0.0009;        4 r0 A& p6 B# e
   
9.     end/ R' h/ s% z# Q4 S- `
   
10.     x = x+0.0011;  O! Y7 L: Y6 N: a\" X/ S
    
11.     y = 1;( \9 T) N0 a, g, i' x- Z0 E
    
12. end: ]; W, p\" x6 M( y. F0 t) u6 T
    
13. s
    6 v$ K% W; X+ }7 a0 ]6 N
14. toc
    + ^5 T# {5 F\" O. }) V2 W
15. 
    , ~! C7 v3 @% L
16. s =
    % e% c# n/ G\" q* m3 `, n\" `. D
17. 1 h4 ?* V7 h3 k8 M5 X+ |
    
18.   1.0086e+006
    6 `* q' M\" {+ s' B! B
19. 
    7 @* A8 J1 ~$ K' N* B: x$ w0 |! d' w
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:5 b+ y$ w) H4 n  p3 y+ M8 t
   
2. t=sys::clock();- Y, e- T+ a\" s$ z4 F1 ]
   
3. s=0,x=0,
   ! T8 b. i; k% r\" J  D' u. H
4. while{x<=1, //while循环算法；
   + s/ ^5 |* |% T\" ~0 r8 F0 ?
5.     y=1, : v3 }! ?$ z% X7 g0 g- t, {
   
6.     while{y<=2, 3 j; m& l; H2 Y3 B: g# I7 S3 y
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), & X7 p* y( Q6 }# X; |: k
   
8.         y=y+0.0009 2 o# t% {2 Y3 Z$ |2 i) x
   
9.     },
   / v( N* Z6 r% H( P* L\" @
10.     x=x+0.0011 . d8 c' D( i' g% q
    
11. },
    : d  E4 E/ n5 b2 O# q  S
12. s;) f8 q\" @, p$ \7 z
    
13. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
6 V4 Q3 e) S) k8 l0.734   //时间，秒

. O/ K/ ^& W1 Z# |我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
8 s' y$ v4 c4 C/ K
2 Q3 m6 x/ i$ X－－－－－－－

2 a  h6 t' O6 l9 b# oForcal中还有一个函数sum专门进行这种计算：

1. mvar:\" |1 e- n+ r6 S: f% r
   
2. t=sys::clock();2 w1 L& e# ~, ?  M/ I
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ) }$ r: H2 H6 A5 V
4. sum["f",0,1,0.0011 : 1,2,0.0009];; \3 u3 f: t1 _! c2 h\" x. W
   
5. [sys::clock()-t]/1000;

复制代码

结果：
: A6 ~1 K8 \# k* A8 b& K  x9 F2 x1008606.64947441
6 t+ G0 P2 |5 @3 H2 k0.719   //时间，秒

sxjm567

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

海水

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

我就在你背后

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

forcal

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

matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   9 @3 b9 Q) b3 O
2. tic
   ( ^) d: Y& U$ j$ J/ A
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   ! _1 L. Z' U  [: N
4. sum(sum(arrayfun(f,x,y)))
   : Q' w6 v' [5 Z, ]( x+ u
5. toc! I2 }( y1 m: Z2 W+ g: A! N
   
6. \" \2 K1 U; F! v7 W& G. b\" m
   
7. ans =$ }5 [3 |! U  G4 \# ]
   
8. 4 y  A+ p# A. L
   
9.   1.0086e+006) }) v' D  |6 f
   
10. 
    * I( _& p% i+ f; G, d, j/ I
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];* |1 g! g4 s9 E: X& Y5 u( L7 a9 W
   
4. mvar:
   
5. ! Z; U- Z/ z& \. j4 K2 z& i! v
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 7 F' P$ x' j% w# o0 ]7 {. n3 @
   
8. t=clock(),- j9 G! P3 i2 X\\" S& h
   
10. oo{
    
11. 8 w/ o% A' k! \\\" p$ M3 z6 B, {
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
13.   \9 P  T6 T\\" T1 F
14.     arrayfun[HFor("f"),x,y].Sum[]
    
15. # ~4 w; K$ _1 b) m! z
16. };
    
17. ! U* B- Z) C% e7 C& C
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">* |1 g! g4 s9 E: X& Y5 u( L7 a9 W</font><br /> mvar:<br /> <span style=\\"display:none\\">! Z; U- Z/ z& \. j4 K2 z& i! v</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">7 F' P$ x' j% w# o0 ]7 {. n3 @</font><br /> t=clock(),<font class=\\"jammer\\">- j9 G! P3 i2 X\\" S& h</font><br /> oo{<br /> <span style=\\"display:none\\">8 w/ o% A' k! \\\" p$ M3 z6 B, {</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">  \9 P  T6 T\\" T1 F</span>     arrayfun[HFor("f"),x,y].Sum[]<br /> <span style=\\"display:none\\"># ~4 w; K$ _1 b) m! z</span> };<br /> <span style=\\"display:none\\">! U* B- Z) C% e7 C& C</span> [clock()-t]/1000;

结果：
1008606.64947441
6 X8 x: R0 Z. c3 F8 r  X0.735  秒

5 M1 s+ C/ d4 E, j, Z可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
! _4 }. J4 S0 |2 e& e* Y; C) W/ q
－－－－－－－－

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