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

forcal

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

1. //用C++代码描述为：
   5 |3 T( z8 h) ?4 B4 h
2. s=0.0;
   ! n5 Y' V6 E/ q4 W' V4 r
3. for(x=0.0;x<=1.0;x=x+0.0011) 4 R! Y0 x7 q! [8 a( n! B
   
4. {, ~( a. K: s5 e\" n' o7 G
   
5.    for(y=1.0;y<=2.0;y=y+0.0009)6 f# b. Z% O' x6 e! r$ I2 S& m& \
   
6.    {8 x: v/ C; }, W
   
7.      s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' P7 ^# W; y. D8 ^6 ]
8.    }
   1 K3 T0 d8 Z7 n7 i2 I
9. }

复制代码

Matlab代码：

1. tic
   ; L; g; m8 i* _) _2 C8 K. {
2. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);, X5 {, l: u) V
   
3. s=sum(sum(cos(1-sin(1.2*x.^(y/2)+cos(1-sin(1.2*y.^(x/2))))))): \2 T! y# f/ r5 E
   
4. toc/ Q7 S; A& w0 ]/ s2 p
   
5. ( z\" `. N$ q1 e7 O
   
6. s =  w2 I+ m2 L! ~) G
   
7. 9 w. `# Y4 T3 C) E( n! O
   
8.   1.0086e+0065 c$ H. c3 o\" w+ V8 `( p% |* K
   
9. 
   6 s+ g7 N! ^8 A
10. Elapsed time is 0.561108 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];# G1 y8 x' H, b5 G: u
   
4. mvar:' [4 Q; C0 m; j
   
6. t=clock(),
   
7. 0 ], g5 B% s5 ^
8. oo{
   
9. % G3 w3 v% [7 T! i0 u
10.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. 3 f9 N+ Z! [( e. R8 a( U
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. ; B9 v/ c* l: r\\" \/ e0 {
14. };/ n3 n\\" c8 v) M
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\"># G1 y8 x' H, b5 G: u</font><br /> mvar:<font class=\\"jammer\\">' [4 Q; C0 m; j</font><br /> t=clock(),<br /> <span style=\\"display:none\\">0 ], g5 B% s5 ^</span> oo{<br /> <span style=\\"display:none\\">% G3 w3 v% [7 T! i0 u</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">3 f9 N+ Z! [( e. R8 a( U</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\\">; B9 v/ c* l: r\\" \/ e0 {</span> };<font class=\\"jammer\\">/ n3 n\\" c8 v) M</font><br /> [clock()-t]/1000;

结果：
! B* v, `$ G6 v3 ?6 m1008606.64947441
0.641
2 Z9 V+ H. E) K. Z
) Y/ q, y+ A/ m+ `, HForcal比Matlab稍慢些。
" ^% C4 |+ t( r4 N8 Z' v4 B
3 W! D5 Z1 x9 n5 M. Q( J* u  j－－－－－－－－－－
% L2 A! E& E# M' s3 D# |
* L* ^/ K: }9 p7 M: A; E" O再看循环效率。
& Q3 M" Y& M6 W( o& P& v
! W3 {; ~$ y  tMatlab代码：

1. tic
     U2 A8 c% d+ H* F
2. s=0;! i; v6 l4 c) K+ Q\" Y( `
   
3. x = 0;0 ?% A2 Q% h$ Y4 p
   
4. y = 1;
   3 f: O. ]# `& W' E
5. while x<1    + Y9 \) u: ^5 D/ o' S8 g  D
   
6.     while y<2;        6 R\" @# `! k' |$ ?) j
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));% ?9 t* f\" s5 Q- G
   
8.         y = y+0.0009;        
   & b! j7 Q8 B# B\" k& v
9.     end
   ; F' O4 d$ |4 `, U, G8 `4 y
10.     x = x+0.0011;
    & a: S2 E- _: [, y& B5 j# D/ y& Z
11.     y = 1;+ d1 E9 F2 O) ?\" B; m7 B\" K
    
12. end3 l* F6 N$ G\" M5 i% P
    
13. s% {( V7 o; M8 l& D$ h9 a\" J
    
14. toc0 G8 g# z: |0 R* h2 W
    
15. 8 u2 y: J+ l+ Q, e* u3 o
    
16. s =
    / e\" K. r8 e9 K$ u' z) Z: x
17. 5 C4 ?  f8 o3 A% `0 ^
    
18.   1.0086e+006
    3 F7 Y. [: \/ I# H% p# V9 U* [
19. 
    8 A/ ]5 O: ^$ D8 t# i
20. Elapsed time is 0.933513 seconds.

复制代码

Forcal代码：

1. mvar:
   ' P, n) i' W# V& r% ]6 [
2. t=sys::clock();' B& y0 L! w. [. p9 Q. M6 g
   
3. s=0,x=0,
   ; Z) N% h  j1 }# Z
4. while{x<=1, //while循环算法； + w* M\" i; S' Z% ^! F% D& X
   
5.     y=1,
   \" S8 t  K& r% X* l3 o
6.     while{y<=2, / ?- j- r' G/ R# m3 T
   
7.         s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ) _9 G# v$ l- r\" i1 N- G
8.         y=y+0.0009 3 I' y: R7 C8 ?: h, u3 c8 T
   
9.     },
   * }+ K+ K! `7 y. d; `$ r
10.     x=x+0.0011
    - L2 [4 C/ n  r
11. },
    & N) w& o* j4 }
12. s;
    2 o, w# g, k) p! {! v: O* l# g
13. [sys::clock()-t]/1000;

复制代码

结果：
' v) R" K, }4 Y. ^! F! j1008606.64947441
3 |* Y( J) N" S5 q0.734   //时间，秒

我很奇怪，在这个例子中，matlab的JIT加速器为什么没有起作用？
+ z3 u) ]  x6 O% b+ Y
/ ?" ]" T. I# ?7 ]" S3 C& D: P－－－－－－－
7 O1 H# s- k) a6 f
: e% ~5 f' X; f. v1 s5 NForcal中还有一个函数sum专门进行这种计算：

1. mvar:% P; s( o: G3 N# ]/ j5 ^
   
2. t=sys::clock();; y. l' c' d) P( ~& c5 _
   
3. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   6 s* p. F( }3 z$ J$ l6 K& P
4. sum["f",0,1,0.0011 : 1,2,0.0009];
   $ Z- i( r0 n7 i% H
5. [sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
9 _# n  \, ?3 S2 q0.719   //时间，秒

forcal

说到代码矢量化，就不能不提到大名鼎鼎的arrayfun函数，现在再来看看它的表现。
5 @* Y9 f0 J# e" o3 `2 r
, s& y. ~9 V" [matlab代码：

1. f=@(x,y)cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ' D# \8 m/ m, ?9 r
2. tic
   6 e' {. j: E  \0 `2 w( a! l
3. [x,y]=meshgrid(0:0.0011:1,1:0.0009:2);7 y* j9 f% i3 _/ m- ^
   
4. sum(sum(arrayfun(f,x,y)))
   # a: Y, u. g  a) R) d& @
5. toc8 R\" {1 `: q& Z; h$ K
   
6. \" x1 U\" H$ R' Z) G7 @* U9 O
   
7. ans =
   7 O& W, ^1 ~) U- w+ D+ {4 z
8. , {# `( j9 A( S/ b4 p
   
9.   1.0086e+006
   & d- t4 ]3 A, D0 L2 z
10. 
    7 N& [! W' K9 P
11. Elapsed time is 7.339923 seconds.

复制代码

Forcal代码：

2. !using["math","sys"];
   
3. 2 s- i9 M* m' W3 n$ Z6 Z
4. mvar:
   
5. * ^9 G& g2 `( n9 l
6. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); * U2 `! `- ]. n6 Y( d( ], f
   
8. t=clock(),
   
9. / a9 q6 q! G% q' J7 h
10. oo{
    
11. / Z\\" f+ w0 G0 t+ \5 r
12.     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],\\" v; |, R3 B4 A% O* w
    
14.     arrayfun[HFor("f"),x,y].Sum[]. g, \, g. d5 f: C$ K' z* M  v
    
16. };! T  ~9 c+ N9 i. I, R1 J
    
18. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">2 s- i9 M* m' W3 n$ Z6 Z</span> mvar:<br /> <span style=\\"display:none\\">* ^9 G& g2 `( n9 l</span> f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); <font class=\\"jammer\\">* U2 `! `- ]. n6 Y( d( ], f</font><br /> t=clock(),<br /> <span style=\\"display:none\\">/ a9 q6 q! G% q' J7 h</span> oo{<br /> <span style=\\"display:none\\">/ Z\\" f+ w0 G0 t+ \5 r</span>     ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">\\" v; |, R3 B4 A% O* w</font><br />     arrayfun[HFor("f"),x,y].Sum[]<font class=\\"jammer\\">. g, \, g. d5 f: C$ K' z* M  v</font><br /> };<font class=\\"jammer\\">! T  ~9 c+ N9 i. I, R1 J</font><br /> [clock()-t]/1000;

结果：
8 F  X3 {/ V2 S$ x, z# h4 K4 I. D1008606.64947441
$ A8 n  s0 i% Q* V0.735  秒

可以看出，Forcal的arrayfun函数远快于matlab的arrayfun函数。
. q0 n4 o1 M, j2 d, p7 A8 `1 I/ x1 C
" J% y# F. n! Q1 q－－－－－－－－

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

我就在你背后

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

海水

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

sxjm567

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