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

forcal

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

1. //用C++代码描述为：
   * D- `! D6 r\" p% j; F
2. s=0.0;
   9 C- p1 _' _  g) x
3. for(x=0.0;x<=1.0;x=x+0.0011)
   # @+ v+ }* H\" A6 L3 O+ O
4. {) A1 n( p- U/ N+ z  o) W
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)9 N% S3 {7 ?: g% @
   
6.    {
   6 G6 Q7 b# A* p# `) l2 ], g
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));& a\" _! ?/ q. M
   
8.    }
     ]- l\" S1 D7 |2 Z5 p4 S& H) h
9. }

复制代码

Matlab代码：

1. tic
   \" h/ P2 T* u2 b7 e\" H
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);/ F+ e  G: V1 z( h, W
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2)))))))
   9 s% j0 I& Q6 r
4. toc& {6 O5 C' Q7 u& R- J
   
5. 4 X% {% \% z$ p- m5 U) y
   
6. s =9 t% q. L% x, C( {! T0 B1 X1 J
   
7. 
   ) b* P) d# c8 u1 W4 z; f  R
8.   1.0086e+006
   5 s. @% ~' ^1 `( s
9. , ]- H! h: h, R3 H' P
   
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];% Y$ U5 D) ~0 Q# L
   
4. mvar:  v\\" T+ u, o0 z: a# h
   
6. t=clock(),
   
7. / [# O. x, t, `+ b* d' M
8. oo{
   
9. - I5 `+ U' v* h3 j5 a- F' q
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. + n- j6 }; W\\" C6 @
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. 2 E1 o7 f! c3 ^8 n
14. };
    
15. + }9 Z) U1 _8 U8 F3 O$ p! o
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">% Y$ U5 D) ~0 Q# L</font><br /> mvar:<font class=\\"jammer\\">  v\\" T+ u, o0 z: a# h</font><br /> t=clock(),<br /> <span style=\\"display:none\\">/ [# O. x, t, `+ b* d' M</span> oo{<br /> <span style=\\"display:none\\">- I5 `+ U' v* h3 j5 a- F' q</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">+ n- j6 }; W\\" C6 @</span>     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\\">2 E1 o7 f! c3 ^8 n</span> };<br /> <span style=\\"display:none\\">+ }9 Z) U1 _8 U8 F3 O$ p! o</span> [clock()-t]/1000;

结果：
% Z0 J4 t  G; [* Y/ Q& v1008606.64947441
0.641
9 Z. z  K7 a* z7 B
  q% k! g- P/ M  ?' U9 t" aForcal比Matlab稍慢些。
( l  |' R6 R* M" c7 _
－－－－－－－－－－

再看循环效率。

Matlab代码：

1. tic! `/ \/ X7 b( u& G
   
2. s=0;
   9 A0 |0 l/ ~: s6 ]; o% @8 Y
3. x = 0;2 Z+ Y+ \! z$ Z\" t3 x
   
4. y = 1;
   2 L2 u: \' j4 }& @\" y\" T5 p* Z9 U/ e
5. while x<1   
   / g' O- a' e+ n  O3 i7 K
6.     while y<2;        
   7 P5 ]2 {% b, e% l
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   # B$ z  i2 T' B  z$ m
8.         y = y+0.0009;        ! F) y3 L. K& Q2 w: p1 K/ y, l
   
9.     end
   1 q  u0 M1 I, j0 j3 R
10.     x = x+0.0011;4 R9 C! Q0 p9 I$ F
    
11.     y = 1;
    ; |' O# E1 k: E0 g) b3 W
12. end
    - ?4 Y! f6 W\" ^% I: K* Z
13. s1 q: n- G' v3 a; r! G3 `  p9 s
    
14. toc7 x* K) d  K6 E  I' ]\" Q! D5 q
    
15. 
    / O\" b4 I/ Y. W) G! F
16. s =6 L9 h/ n! f' |0 _
    
17. 
    . q\" G8 _* j6 X# X- ^! H. I
18.   1.0086e+006
    - c5 ?. ~# J# I
19. 2 k  d0 l# p: r/ K* Y! d
    
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   5 h* F/ k% _1 H
2. t=sys::clock();2 Y) s! [1 I# n% G: O9 _) B: a
   
3. s=0,x=0,
   4 S7 ]; T. h\" l4 _\" Y9 @8 v
4. while{x<=1, //while循环算法；
   8 X* [! f% J7 [& w) |9 N
5.     y=1, 2 y3 d+ I0 z6 X* P
   
6.     while{y<=2,
   1 j- ^: t. U# B2 N5 s- O
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   1 ]% I9 C\" o: V, v* K; T
8.         y=y+0.0009 % w2 [/ Y6 ]; ^
   
9.     }, 6 P5 ?\" i9 p. A( Z* A# G1 s\" b
   
10.     x=x+0.0011 ' A* q5 |* b6 z4 t7 f& i9 w+ u
    
11. }, * X* t- L6 o; _1 e9 g; L
    
12. s;
    - W0 a% h3 ~2 J& E; Q) c2 N8 H
13. [sys::clock()-t]/1000;

复制代码

结果：
+ Z+ ^) Q/ n: S9 I& b1008606.64947441
0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
0 L) R$ ~9 h) Z
－－－－－－－
, m4 l  D! O; o! P: ^- i
Forcal中还有一个函数sum专门进行这种计算：

1. mvar:
     k0 S( h5 N. ?  J4 [
2. t=sys::clock();/ u) t  d  c8 ~, ?\" v/ E0 m
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   1 E, ^9 ~2 h) L6 p& W
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   # U/ s/ E0 C( ^0 a  @9 y
5. [sys::clock()-t]/1000;

复制代码

结果：
8 A0 a4 r0 o2 Y/ o5 R1008606.64947441
) ]+ d4 c0 I2 R6 a0.719   //时间，秒

forcal

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

" E6 R  W9 V, t. B& ]) `7 p$ T: C, Cmatlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ( e( H# d. h& h% t, L
2. tic. X5 I( s# [3 v1 r  A$ ~
   
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);
   & d$ U5 b; k- t9 H5 m
4. sum(sum(arrayfun(f,x,y)))$ _3 F: e: y: P/ k9 M; F
   
5. toc% z4 q\" j  u6 F2 A! u% @, K2 j
   
6. 0 ^: S  ?( a* h8 ?. s
   
7. ans =! p2 H* T* ]9 K/ b* _; H- L
   
8. ! y8 m/ Y4 c- H4 N- o1 c
   
9.   1.0086e+006
   & d& q) u& P$ w  j
10. 
    & c( @& |, S5 J) T7 T$ ~
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];# j6 [0 b  m9 m2 X
   
4. mvar:6 C$ [4 U& B# S, l9 {- z
   
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); ) ]* C) N4 n0 t, w
   
8. t=clock(),
   
9. ) f' _3 ]+ D  x# |* X2 e
10. oo{1 s6 |/ C$ Z6 ?
    
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],5 ^  F4 o3 O+ l. ]( N- c\\" X5 ]
    
14.     arrayfun[HFor("f"),x,y].Sum[]2 f2 C: ^# t% B8 r8 E+ Z+ e
    
16. };. W$ Q1 m  f5 `' ~) \+ D
    
18. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\"># j6 [0 b  m9 m2 X</font><br /> mvar:<font class=\\"jammer\\">6 C$ [4 U& B# S, l9 {- z</font><br /> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">) ]* C) N4 n0 t, w</font><br /> t=clock(),<br /> <span style=\\"display:none\\">) f' _3 ]+ D  x# |* X2 e</span> oo{<font class=\\"jammer\\">1 s6 |/ C$ Z6 ?</font><br />     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">5 ^  F4 o3 O+ l. ]( N- c\\" X5 ]</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">2 f2 C: ^# t% B8 r8 E+ Z+ e</font><br /> };<font class=\\"jammer\\">. W$ Q1 m  f5 `' ~) \+ D</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.735  秒
, k$ [3 A: S0 m0 Y+ i
可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。

－－－－－－－－
% U( n1 z7 [, Y, |7 |% L% b
" L- V: b' W) `! B; R7 @5 G从这里似乎可以看出，matlab若借助于arrayfun函数计算三重及以上积分，其效率将远远落后于Forcal。当然，Forcal不使用arrayfun函数计算三重及以上积分。

我就在你背后

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

海水

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

sxjm567

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