Forcal数学库FcMath：以矩阵运算为基础

forcal

    FcMath32W.dll是一个Forcal数值计算扩展动态库，该库以线性代数特别是矩阵运算为基础。

    在FcMath中的函数是通过二级函数命名空间“math”输出的，所有函数均具有类似“math::array(...)”的格式，都是实数函数。使用!using("math");可简化FcMath中的函数访问。

    FcMath32W.dll需要FcData32W.dll的支持。FcData32W.dll要先于FcMath32W.dll加载。

    FcMath库的数组是C格式的，元素序号是基于0的。可以使用函数sys::rearray在Forcal数组（C数组格式）和Fortran数组之间进行转换。

    一般，若FcMath函数返回一个对象，则在oo函数中将返回临时对象，否则返回一般对象；临时对象由oo函数进行管理，一般对象须用函数delete销毁。故若没有特殊的原因，建议在oo函数中使用FcMath函数！若一般对象没有及时用delete销毁，则其将常驻内存，消耗内存资源；可用FcData的函数DelAllFCD()销毁所有对象，释放内存资源，或者在程序退出时自动销毁所有对象。

    FcMath库函数具有内存消耗低、执行效率高、代码简洁、实用性强的特点。

    FcMath库中所用的算法或许不是最好的，如果您有好的算法，可以方便地进行替换，提升FcMath的性能。

    FcMath库可用于开发极致性能的应用程序，是熟悉C/C++、Fortran的数学爱好者的极佳的练手工具，同时也期望对一般的数值计算用户提供越来越多的方便。

    限于作者水平，期待与朋友们共同完善FcMath！如果您有什么好的算法，任何改进的意见或建议，请与作者联系。

& g- w5 J; Y2 m: m

    详细说明：Forcal数值计算扩展动态库FcMath

    源代码下载：http://www.forcal.net/xiazai/forcal9/forcal9code.rar

forcal

例子1代码：
# y! }/ j+ i+ D' [* N

2. !using["math"];
   
3. * P/ n, T0 J) X/ y, E
4. mvar:
   
5. ) z3 h, E7 q' ?5 V! j# O
6. oo{                      //一般在oo函数中调用FcMath函数
   
7. 2 A( r1 j8 ]& n; g
8.   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化
   
9. / i9 B8 a1 k8 c* I; B, f1 U
10.   a.outm(),              //输出矩阵a$ k& ^0 j$ _: l) t; a
    
12.   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出' i2 G) M  U8 }
    
14.   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出
    
15. . ?5 U5 N) S/ Q! H8 x: i+ C% z
16.   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出% x, |# ]7 `) w( p  ^, G5 G
    
18. };
    
19.   n- ~2 u7 Y6 P4 g7 ~: X& e$ ]

!using["math"];<br /> <span style=\\"display:none\\">* P/ n, T0 J) X/ y, E</span> mvar:<br /> <span style=\\"display:none\\">) z3 h, E7 q' ?5 V! j# O</span> oo{                      //一般在oo函数中调用FcMath函数<br /> <span style=\\"display:none\\">2 A( r1 j8 ]& n; g</span>   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化<br /> <span style=\\"display:none\\">/ i9 B8 a1 k8 c* I; B, f1 U</span>   a.outm(),              //输出矩阵a<font class=\\"jammer\\">$ k& ^0 j$ _: l) t; a</font><br />   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">' i2 G) M  U8 }</font><br />   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出<br /> <span style=\\"display:none\\">. ?5 U5 N) S/ Q! H8 x: i+ C% z</span>   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">% x, |# ]7 `) w( p  ^, G5 G</font><br /> };<br /> <span style=\\"display:none\\">  n- ~2 u7 Y6 P4 g7 ~: X& e$ ]</span>

结果：

1.        0.211319   4.91638e-002       0.144638       0.153259       0.8526159 ~' u0 c5 }- ~4 B9 y5 x2 p, `
   
2.        0.630646       0.927048       0.440308       0.162857       0.556854
   2 ~7 J3 }7 ]+ y, i; ~  |) f/ N
3.         0.43309        0.34552       0.563919       0.937164       0.209641
   9 y% k: I* Z7 ^: K/ p. w6 b
4.        0.603271       0.727676       0.130951   5.35736e-002       0.1979375 E/ @: ~$ e( ^9 S0 Q: x
   
5.        0.576004       0.747589   1.17645e-002       0.363892       0.2807774 g5 P% {! v8 X\" N! I# k
   
6.        0.646454       0.381088        0.58551        0.26387        0.93692
   , G$ M  u) r3 c) @3 q3 l\" m
7. 
   9 m6 f. D4 F. R+ s3 k6 }- h! B\" }
8.        0.153259
   ! i( k1 U9 d. |' ?2 O3 M0 x# U
9.        0.162857
   9 z  x& h8 h7 f% ]( ?- u' }
10.        0.937164
    % p' Y' w: t\" T7 j# S& k
11.    5.35736e-002
    4 `5 l9 e+ l; I0 e
12.        0.363892
    / A* p$ X  P, l$ e) R
13.         0.26387
    : [# P: S, |$ x* Q\" \9 E, l) C\" H
14. 8 r1 [/ K# K( u* V1 {3 P- r
    
15.        0.603271       0.727676       0.130951   5.35736e-002       0.197937# |4 e* U. |$ J* u9 ~  n4 L
    
16. $ Y  C: E- [1 g4 |) s2 d
    
17.        0.130951   5.35736e-002+ n2 N' d1 {+ D' k% y, X/ d
    
18.    1.17645e-002       0.363892. ]1 o, y) r7 H; s8 H0 I4 ?0 q% p
    
19.         0.58551        0.26387$ V* w4 }0 z( M2 o3 Q
    
20. / X$ A6 G& ]( A( D
    

复制代码

" {1 v# t' V6 A/ O6 _例子2代码：

2 y# b" g+ Y* l  `9 `  B

2. f(x1,x2,x3,y1,y2,y3)=      //函数定义
   
3. 6 w0 [$ E1 z3 r% p
4. {
   
5. \\" q( W' w\\" u9 I% A  f
6.     y1=x1*x1+x2*x2+x3*x3-1.0,
   
7. 5 M8 i4 h8 V! V+ y
8.     y2=2.0*x1*x1+x2*x2-4.0*x3,; o\\" O$ |4 k4 B: B# r
   
10.     y3=3.0*x1*x1-4.0*x2+x3*x3: E! q6 L0 j$ G% F* ]9 B6 M
    
12. };
    
13. , i9 G! E3 G7 n) ~4 ?/ F( W& E2 z
14. !using["math","sys"];
    
15. % ?$ X- L# K. I! J' N! u8 R
16. mvar:
    
17. \\" E& n/ B  Y1 a+ [+ s2 w/ s5 k\\" g
18. oo{
    
19. % W\\" f3 }4 x3 r
20.   x=array(3),- |- a& h2 N/ F\\" E# ?
    
22.   x.SA[0 : 1,1,1],       //设置初值为1,1,1
    
23.   ^\\" l1 I\\" m1 }
24.   i=netn[HFor("f"),x],   //拟牛顿法解方程- D) N& i. e; F0 {5 D$ j
    
26.   x.outm(),              //输出结果) e. o1 l) `* R' ]% l$ g0 B4 B3 H. p/ m# x
    
28.   i                      //返回迭代次数! p8 T) C' {) d8 o. E3 X
    
30. };
    
31. 0 A  t- L; ~9 w# x) s8 d

f(x1,x2,x3,y1,y2,y3)=      //函数定义<br /> <span style=\\"display:none\\">6 w0 [$ E1 z3 r% p</span> {<br /> <span style=\\"display:none\\">\\" q( W' w\\" u9 I% A  f</span>     y1=x1*x1+x2*x2+x3*x3-1.0,<br /> <span style=\\"display:none\\">5 M8 i4 h8 V! V+ y</span>     y2=2.0*x1*x1+x2*x2-4.0*x3,<font class=\\"jammer\\">; o\\" O$ |4 k4 B: B# r</font><br />     y3=3.0*x1*x1-4.0*x2+x3*x3<font class=\\"jammer\\">: E! q6 L0 j$ G% F* ]9 B6 M</font><br /> };<br /> <span style=\\"display:none\\">, i9 G! E3 G7 n) ~4 ?/ F( W& E2 z</span> !using["math","sys"];<br /> <span style=\\"display:none\\">% ?$ X- L# K. I! J' N! u8 R</span> mvar:<br /> <span style=\\"display:none\\">\\" E& n/ B  Y1 a+ [+ s2 w/ s5 k\\" g</span> oo{<br /> <span style=\\"display:none\\">% W\\" f3 }4 x3 r</span>   x=array(3),<font class=\\"jammer\\">- |- a& h2 N/ F\\" E# ?</font><br />   x.SA[0 : 1,1,1],       //设置初值为1,1,1<br /> <span style=\\"display:none\\">  ^\\" l1 I\\" m1 }</span>   i=netn[HFor("f"),x],   //拟牛顿法解方程<font class=\\"jammer\\">- D) N& i. e; F0 {5 D$ j</font><br />   x.outm(),              //输出结果<font class=\\"jammer\\">) e. o1 l) `* R' ]% l$ g0 B4 B3 H. p/ m# x</font><br />   i                      //返回迭代次数<font class=\\"jammer\\">! p8 T) C' {) d8 o. E3 X</font><br /> };<br /> <span style=\\"display:none\\">0 A  t- L; ~9 w# x) s8 d</span>

; y1 @. I! e/ U* B结果：
  0.785197       0.496611       0.369923
, E' k3 ]# G1 k. }9 b' ]% c8 ?

forcal

效率测试：
+ K3 q+ ^3 C& B! Y, l; I0 Ssimwe的网友lin2009 的matlab代码：

1. clear all
   # a5 n6 }5 @# G4 B$ W6 D( ~/ x8 n' L
2. clc; n$ P) K: O7 s, U
   
3. tic
   7 _- s( e7 w. g$ s/ O9 y$ {& \
4. k = zeros(5,5); % //生成5×5全0矩阵8 }/ G0 ]5 {) P' n
   
5. % 循环计算以下程序段100000次：4 t1 c, U( L2 x+ j* c4 a4 Z+ _& f
   
6. for m = 1:100000\" g) m( `1 q$ M, n% N7 ~
   
7.     a = rand(5,7);
   \" [' r- h+ W0 J! V7 W
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   * j' j: G, @# Y
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);6 N3 e! d8 w& t( ^
   
10. end
    : j' l: E* w4 c
11. k, v! O/ H& c  e
    
12. toc
    5 h, ?' S8 W  j, F8 y9 F4 m

复制代码

Forcal代码：
' D9 e+ \1 h" s( P7 u; y" d
运行稍快的代码，比matlab约快10％吧？

& G5 k$ u+ X3 p0 ~  B

2. !using["math","sys"];
   
3. 9 A7 I- h$ \; T' v0 g+ n
4. mvar:
   
5. 1 q) ~8 u  a; u+ I! g9 u
6. t0=clock(),) H3 v# l! k8 `) X! `% i7 s
   
8. oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为07 u: j8 E3 x( J  Y. q- \' Y
   
10. i=0,(i<1000 00).while{ //循环计算1000 00次
    
11. / N5 Z: Q4 c* H
12.   oo{
    
13. 4 b5 k+ I\\" c3 ~3 r) V* C, h0 \
14.     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
    
15.   B9 {% d% j/ e$ b# E
16.     k.oset[k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)] //计算k=k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)
    
17. / X  H% ]& Z; r4 }; u6 l' v
18.   },4 p+ C5 x: O1 \) H9 \
    
20.   i++
    
21. 0 ]  v) F% C; X* B5 H
22. },/ \9 \8 ]# O- G3 c/ |) _
    
24. k.outm(),             //输出矩阵k，然后销毁k
    
25. \\" r: Q+ W5 v  [0 a& f+ m  b\\" A
26. [clock()-t0]/1000;    //得到计算时间，秒

!using["math","sys"];<br /> <span style=\\"display:none\\">9 A7 I- h$ \; T' v0 g+ n</span> mvar:<br /> <span style=\\"display:none\\">1 q) ~8 u  a; u+ I! g9 u</span> t0=clock(),<font class=\\"jammer\\">) H3 v# l! k8 `) X! `% i7 s</font><br /> oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0<font class=\\"jammer\\">7 u: j8 E3 x( J  Y. q- \' Y</font><br /> i=0,(i<1000 00).while{ //循环计算1000 00次<br /> <span style=\\"display:none\\">/ N5 Z: Q4 c* H</span>   oo{<br /> <span style=\\"display:none\\">4 b5 k+ I\\" c3 ~3 r) V* C, h0 \</span>     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化<br /> <span style=\\"display:none\\">  B9 {% d% j/ e$ b# E</span>     k.oset[k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)] //计算k=k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)<br /> <span style=\\"display:none\\">/ X  H% ]& Z; r4 }; u6 l' v</span>   },<font class=\\"jammer\\">4 p+ C5 x: O1 \) H9 \</font><br />   i++<br /> <span style=\\"display:none\\">0 ]  v) F% C; X* B5 H</span> },<font class=\\"jammer\\">/ \9 \8 ]# O- G3 c/ |) _</font><br /> k.outm(),             //输出矩阵k，然后销毁k<br /> <span style=\\"display:none\\">\\" r: Q+ W5 v  [0 a& f+ m  b\\" A</span> [clock()-t0]/1000;    //得到计算时间，秒

在我的电脑上运行时间为3.344秒。

; T" l* c( ]) V+ j" A比较好看些的代码，似乎也比matlab稍快吧？
& [  d( i, q  m* D

2. !using["math","sys"];
   
3. * v8 p$ m1 a: o0 E; ^* V
4. (:t0,k,i,a,b)=5 o8 N- w8 v% i/ Y* q
   
6. {. ^$ O3 ~  p8 r  O9 N  U' {
   
8.   t0=clock(),; D. V6 h; e\\" p6 M- X. v6 l
   
10.   k=zeros[5,5],
    
11. % \, c1 |' n4 b( H  P
12.   i=0,(i<1000 00).while{& i% c  K! X, o2 P* F9 _; c
    
14.     oo{
    
15. / K! m( b\\" G1 J7 S\\" y
16.       a=rand[5,7], b=rand[7,5],
    
17. \\" k1 `( }$ b# X0 @' l' ]3 A
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
    
19. $ }) v( \3 Y8 w/ m0 v! n
20.     },
    
21. 2 u3 w* X) q8 r- h7 l3 F
22.     i++5 ]# U& X6 s' J. z3 G) l* ?
    
24.   },% t& y+ R* H$ V: A( P9 L1 x  n
    
26.   k.outm().delete(),0 l- ^6 i9 B- e+ q  J5 w' X9 B
    
28.   [clock()-t0]/10004 z4 K6 D( y4 h2 Q4 o6 ]
    
30. };

!using["math","sys"];<br /> <span style=\\"display:none\\">* v8 p$ m1 a: o0 E; ^* V</span> (:t0,k,i,a,b)=<font class=\\"jammer\\">5 o8 N- w8 v% i/ Y* q</font><br /> {<font class=\\"jammer\\">. ^$ O3 ~  p8 r  O9 N  U' {</font><br />   t0=clock(),<font class=\\"jammer\\">; D. V6 h; e\\" p6 M- X. v6 l</font><br />   k=zeros[5,5],<br /> <span style=\\"display:none\\">% \, c1 |' n4 b( H  P</span>   i=0,(i<1000 00).while{<font class=\\"jammer\\">& i% c  K! X, o2 P* F9 _; c</font><br />     oo{<br /> <span style=\\"display:none\\">/ K! m( b\\" G1 J7 S\\" y</span>       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\">\\" k1 `( }$ b# X0 @' l' ]3 A</span>       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<br /> <span style=\\"display:none\\">$ }) v( \3 Y8 w/ m0 v! n</span>     },<br /> <span style=\\"display:none\\">2 u3 w* X) q8 r- h7 l3 F</span>     i++<font class=\\"jammer\\">5 ]# U& X6 s' J. z3 G) l* ?</font><br />   },<font class=\\"jammer\\">% t& y+ R* H$ V: A( P9 L1 x  n</font><br />   k.outm().delete(),<font class=\\"jammer\\">0 l- ^6 i9 B- e+ q  J5 w' X9 B</font><br />   [clock()-t0]/1000<font class=\\"jammer\\">4 z4 K6 D( y4 h2 Q4 o6 ]</font><br /> };

在我的电脑上运行时间为3.579秒。
# x6 a3 X% _" Y( E; E% a
该例子的理论结果是每个元素均为275000。
+ _! e- `) D6 T9 W: j9 p% F& q
' q; O7 E/ i3 {$ ^% ^9 a& Y我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。
1 q; s9 J; M$ ]

forcal

继续例子，大家看有什么问题吗？

2. !using["math"];
   
3. 6 R2 p$ e# k! p
4. mvar:+ s/ Y% K1 K; }. L% N1 F' H
   
6. oo{
   
7. ) I0 L. Z; g' u' N* @# C
8.   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],+ r/ d9 j- |/ p; Z
   
10.   a1.outm[5,1,1],/ b0 m- ?$ U! w+ ]4 o& W0 T5 {
    
12.   a2.outm[5,1,1],4 u# V4 \# `4 o& D: ~5 G7 n3 I! X- v
    
14.   a3.outm[5,1,1],
    
15. 8 x6 C/ g) |7 y\\" U) E  p% a& f
16.   a4.outm[5,1,1],
    
17. : a# z0 n0 S+ P- e. c( x% A0 _
18.   a=a1+a2+a3+a4,
    
19. & `, u. ^% F5 `% n
20.   a.outm[5,1,1],
    
21. ) k8 x; h/ N, i& `: N
22.   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]
    
23. $ B# ?; b- R9 c& T; g( E
24. };
    
25. 7 W5 v# F1 j. j& A

!using["math"];<br /> <span style=\\"display:none\\">6 R2 p$ e# k! p</span> mvar:<font class=\\"jammer\\">+ s/ Y% K1 K; }. L% N1 F' H</font><br /> oo{<br /> <span style=\\"display:none\\">) I0 L. Z; g' u' N* @# C</span>   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],<font class=\\"jammer\\">+ r/ d9 j- |/ p; Z</font><br />   a1.outm[5,1,1],<font class=\\"jammer\\">/ b0 m- ?$ U! w+ ]4 o& W0 T5 {</font><br />   a2.outm[5,1,1],<font class=\\"jammer\\">4 u# V4 \# `4 o& D: ~5 G7 n3 I! X- v</font><br />   a3.outm[5,1,1],<br /> <span style=\\"display:none\\">8 x6 C/ g) |7 y\\" U) E  p% a& f</span>   a4.outm[5,1,1],<br /> <span style=\\"display:none\\">: a# z0 n0 S+ P- e. c( x% A0 _</span>   a=a1+a2+a3+a4,<br /> <span style=\\"display:none\\">& `, u. ^% F5 `% n</span>   a.outm[5,1,1],<br /> <span style=\\"display:none\\">) k8 x; h/ N, i& `: N</span>   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]<br /> <span style=\\"display:none\\">$ B# ?; b- R9 c& T; g( E</span> };<br /> <span style=\\"display:none\\">7 W5 v# F1 j. j& A</span>

' N0 ?8 ~2 C, S; o; c说明：
3 I. V4 P$ s) c" w' glinspace(8,12,5)：生成一维数组，共5个元素8～12
a1.outm[5,1,1]：输出**数组a1，连下标一起输出
3 D' q: }6 H2 iSum[a]：设有m维数组（含矩阵）：a(n1,n2,... ...,nm)，则Sum(a,i)对第i维求和，返回一个m-1维数组。若a是一维数组、1×k矩阵、k×1矩阵、i<1或者i>m，则Sum函数返回所有数组元素的和。Sum(a)相当于Sum(a,m)。
1 ]+ P5 }) n* `: W6 ~" ?& k0 e) }
结果（最终求和结果是1320）：
% p5 D7 ^) Z: R
& [5 a& ^' V% @! q(0,0,0,*)              1.            1.            1.            1.            1.
7 P" C! |& w8 t; ]; J(0,0,1,*)              1.            1.            1.            1.            1.
+ C" a4 k/ \% N& c8 _* Q(0,1,0,*)              1.            1.            1.            1.            1.
(0,1,1,*)              1.            1.            1.            1.            1.
" I) O% k) x3 ^1 O5 X, I) K(0,2,0,*)              1.            1.            1.            1.            1.
& G6 V1 d0 ?; w# z; M& _(0,2,1,*)              1.            1.            1.            1.            1.
3 @' b7 \" {# E7 M$ {(1,0,0,*)              2.            2.            2.            2.            2.
(1,0,1,*)              2.            2.            2.            2.            2.
(1,1,0,*)              2.            2.            2.            2.            2.
(1,1,1,*)              2.            2.            2.            2.            2.
(1,2,0,*)              2.            2.            2.            2.            2.
(1,2,1,*)              2.            2.            2.            2.            2.
% ?! V5 @7 ]* U
(0,0,0,*)              3.            3.            3.            3.            3.
7 ]) h9 a  G2 e4 S- X- D, m) [(0,0,1,*)              3.            3.            3.            3.            3.
(0,1,0,*)              4.            4.            4.            4.            4.
(0,1,1,*)              4.            4.            4.            4.            4.
(0,2,0,*)              5.            5.            5.            5.            5.
8 s5 }0 \4 T8 h$ k(0,2,1,*)              5.            5.            5.            5.            5.
# M; B. d2 g- Z0 @1 N! K(1,0,0,*)              3.            3.            3.            3.            3.
(1,0,1,*)              3.            3.            3.            3.            3.
# P6 l3 T3 d2 A1 b4 [( F(1,1,0,*)              4.            4.            4.            4.            4.
+ M. N) p, O' J3 K5 c6 @: r(1,1,1,*)              4.            4.            4.            4.            4.
9 y" ]4 u, h: Q! U+ O(1,2,0,*)              5.            5.            5.            5.            5.
(1,2,1,*)              5.            5.            5.            5.            5.
3 e- A$ {& K0 [( B0 Q+ u- }9 C
6 G5 C6 K; E  U5 @(0,0,0,*)              6.            6.            6.            6.            6.
(0,0,1,*)              7.            7.            7.            7.            7.
(0,1,0,*)              6.            6.            6.            6.            6.
(0,1,1,*)              7.            7.            7.            7.            7.
: {% B$ ]) M5 \5 a& U+ o! l(0,2,0,*)              6.            6.            6.            6.            6.
(0,2,1,*)              7.            7.            7.            7.            7.
' O7 ~4 j1 N- }7 O: J7 F(1,0,0,*)              6.            6.            6.            6.            6.
(1,0,1,*)              7.            7.            7.            7.            7.
/ F* W* X+ _! E- K+ L# H0 Q$ p(1,1,0,*)              6.            6.            6.            6.            6.
(1,1,1,*)              7.            7.            7.            7.            7.
(1,2,0,*)              6.            6.            6.            6.            6.
(1,2,1,*)              7.            7.            7.            7.            7.

(0,0,0,*)              8.            9.           10.           11.           12.
(0,0,1,*)              8.            9.           10.           11.           12.
9 W/ e3 G, T+ ]7 x; E' ~(0,1,0,*)              8.            9.           10.           11.           12.
(0,1,1,*)              8.            9.           10.           11.           12.
(0,2,0,*)              8.            9.           10.           11.           12.
(0,2,1,*)              8.            9.           10.           11.           12.
) y% d- k& G9 `(1,0,0,*)              8.            9.           10.           11.           12.
0 p1 c6 a& ?" \6 x(1,0,1,*)              8.            9.           10.           11.           12.
5 _8 }( [" c7 M3 U$ I. i(1,1,0,*)              8.            9.           10.           11.           12.
(1,1,1,*)              8.            9.           10.           11.           12.
+ j* I4 X; S$ F7 @! T3 [& l(1,2,0,*)              8.            9.           10.           11.           12.
& p) y2 l, C5 |! h1 o1 b" g(1,2,1,*)              8.            9.           10.           11.           12.

4 V0 g; K4 U0 g  Q(0,0,0,*)             18.           19.           20.           21.           22.
(0,0,1,*)             19.           20.           21.           22.           23.
- [: H# B8 K- X(0,1,0,*)             19.           20.           21.           22.           23.
% q$ g8 N# y, M- s) Y" v- u% N3 |(0,1,1,*)             20.           21.           22.           23.           24.
7 j+ D! L, Z% t( U(0,2,0,*)             20.           21.           22.           23.           24.
, R/ H& o6 ]! E5 b/ b4 i1 d- O(0,2,1,*)             21.           22.           23.           24.           25.
3 {/ j+ J# ]6 ]9 l0 F: T5 F(1,0,0,*)             19.           20.           21.           22.           23.
2 p/ R; e& d0 @  Y(1,0,1,*)             20.           21.           22.           23.           24.
6 r, o+ _9 t: P: \% X(1,1,0,*)             20.           21.           22.           23.           24.
(1,1,1,*)             21.           22.           23.           24.           25.
(1,2,0,*)             21.           22.           23.           24.           25.
$ ~$ o1 C  K$ {4 {(1,2,1,*)             22.           23.           24.           25.           26.

) Z$ Y/ E. t) e4 S(0,0,*)            100.          105.
* d  [! m, B9 c9 N(0,1,*)            105.          110.
(0,2,*)            110.          115.
% R9 ~3 e8 P$ ?. J1 x1 a) N(1,0,*)            105.          110.
(1,1,*)            110.          115.
(1,2,*)            115.          120.
( k" R/ y  L4 O8 s
(0,*)            205.          215.          225.
2 |& ?7 a/ \) R) W6 m5 i(1,*)            215.          225.          235.

(0,*)            645.
& f# d7 ^1 [9 b  f) k# A2 X(1,*)            675.
3 ?; u# H5 P. _  m, n7 F
! k( Z# O8 _2 I" ~4 ^. S/ X1320.
4 O  ?# Z' {; f0 i
9 m# Q5 [. {1 l# n

forcal

在matlab中，纯for循环速度最慢。而一半for循环+一半向量化的速度最快，Forcal中也是如此：
& P. B4 u4 m& _  w  `, h7 c

2. !using["math","sys"];
   
3. 8 a  U\\" S  [5 v
4. mvar:
   
5. ' }: K\\" r: _5 ]
6. (:p1,p2,p3,a,b)=- {( o( A. ]% H( G0 b: i
   
8. {2 u, Z% d0 {\\" ^* U7 N% x2 u
   
10.   oo{
    
11. ) I8 V8 W$ p$ I! k\\" L$ f
12.     a=array[1000].rand(),
    
13. 1 {) U( A: |- B& g1 s
14.     b=array[1000].rand(),0 L# ~( ?  H1 y+ t- Z
    
16.     p1=array[1000,1000],
    
17. + j% X1 A* Q5 V, q1 l5 s7 t4 K
18.     p2=array[1000,1000],% F5 s' @2 `. Y( a/ `
    
20.     p3=array[1000,1000],
    
21. ! e: T' |$ h; J5 N& ?- p: p
22.     t0=clock(),& ]7 [2 K/ s\\" K8 N8 z  m/ P
    
24.     ndgrid(a,b,&A,&B),
    
25. ' i6 v9 |! r' E% |7 P/ E' G$ j
26.     p1.=A+B
    
27. % _9 z; Z) W* U1 C
28.   },
    
29. , B: E/ W\\" B$ c! F3 H
30.   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},/ N, k( K2 D3 T7 \) ^
    
32.   lena=FCDLen(a),
    
33. 0 K2 ?+ T/ g$ Y) i
34.   lenb=FCDLen(b),! `* h1 [' r* T+ S* S  f! N+ b
    
36.   t0=clock(),
    
37. - ~3 ^; C: \* G1 e4 E1 j; b% V* G
38.   m = lenb-1, (m>=0).while{  z6 |# X9 Z5 l& `7 I: s5 s9 A+ E! v
    
40.     oo{p2(m,neg) = a+rn[b(m)]},. f) `4 @& e, q3 k( [& c- k- b
    
42.     m--
    
43. 8 T2 ]0 `5 Q9 t5 E$ k
44.   },
    
45. 1 {3 P6 q. ?4 @9 h8 P* p
46.   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},
    
47. 4 S/ G% [# h6 [, c
48.   t0=clock(),1 T. |1 J\\" u4 p
    
50.   m = lenb-1, (m>=0).while{
    
51. $ m$ }$ o3 X2 b4 r3 m
52.     n = lena-1, (n>=0).while{. |4 P9 k0 g3 K% `% R2 o1 `5 p9 v
    
54.       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些* @. E* [4 V/ K. t8 Y
    
56.       A(p3,m,n) = A(a,n)+A(b,m),4 B) n2 d4 O& r* h- v
    
58.       n--
    
59. ; k4 l6 L\\" `' h# v/ c
60.     },
    
61.   m- }3 s* f$ d5 o/ k# i  Y( @+ O
62.     m--: K. `; r( F' ]/ A3 W- _) j& _\\" w\\" ^
    
64.   },
    
65. : W( ^) @3 n6 z% b9 i6 x
66.   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
    
67. ( ?: e\\" O* Z* d9 O
68. };
    
69. 0 l& y\\" S: _6 y) S) n/ n4 g

!using["math","sys"];<br /> <span style=\\"display:none\\">8 a  U\\" S  [5 v</span> mvar:<br /> <span style=\\"display:none\\">' }: K\\" r: _5 ]</span> (:p1,p2,p3,a,b)=<font class=\\"jammer\\">- {( o( A. ]% H( G0 b: i</font><br /> {<font class=\\"jammer\\">2 u, Z% d0 {\\" ^* U7 N% x2 u</font><br />   oo{<br /> <span style=\\"display:none\\">) I8 V8 W$ p$ I! k\\" L$ f</span>     a=array[1000].rand(),<br /> <span style=\\"display:none\\">1 {) U( A: |- B& g1 s</span>     b=array[1000].rand(),<font class=\\"jammer\\">0 L# ~( ?  H1 y+ t- Z</font><br />     p1=array[1000,1000],<br /> <span style=\\"display:none\\">+ j% X1 A* Q5 V, q1 l5 s7 t4 K</span>     p2=array[1000,1000],<font class=\\"jammer\\">% F5 s' @2 `. Y( a/ `</font><br />     p3=array[1000,1000],<br /> <span style=\\"display:none\\">! e: T' |$ h; J5 N& ?- p: p</span>     t0=clock(),<font class=\\"jammer\\">& ]7 [2 K/ s\\" K8 N8 z  m/ P</font><br />     ndgrid(a,b,&A,&B),<br /> <span style=\\"display:none\\">' i6 v9 |! r' E% |7 P/ E' G$ j</span>     p1.=A+B<br /> <span style=\\"display:none\\">% _9 z; Z) W* U1 C</span>   },<br /> <span style=\\"display:none\\">, B: E/ W\\" B$ c! F3 H</span>   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},<font class=\\"jammer\\">/ N, k( K2 D3 T7 \) ^</font><br />   lena=FCDLen(a),<br /> <span style=\\"display:none\\">0 K2 ?+ T/ g$ Y) i</span>   lenb=FCDLen(b),<font class=\\"jammer\\">! `* h1 [' r* T+ S* S  f! N+ b</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">- ~3 ^; C: \* G1 e4 E1 j; b% V* G</span>   m = lenb-1, (m>=0).while{<font class=\\"jammer\\">  z6 |# X9 Z5 l& `7 I: s5 s9 A+ E! v</font><br />     oo{p2(m,neg) = a+rn[b(m)]},<font class=\\"jammer\\">. f) `4 @& e, q3 k( [& c- k- b</font><br />     m--<br /> <span style=\\"display:none\\">8 T2 ]0 `5 Q9 t5 E$ k</span>   },<br /> <span style=\\"display:none\\">1 {3 P6 q. ?4 @9 h8 P* p</span>   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},<br /> <span style=\\"display:none\\">4 S/ G% [# h6 [, c</span>   t0=clock(),<font class=\\"jammer\\">1 T. |1 J\\" u4 p</font><br />   m = lenb-1, (m>=0).while{<br /> <span style=\\"display:none\\">$ m$ }$ o3 X2 b4 r3 m</span>     n = lena-1, (n>=0).while{<font class=\\"jammer\\">. |4 P9 k0 g3 K% `% R2 o1 `5 p9 v</font><br />       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些<font class=\\"jammer\\">* @. E* [4 V/ K. t8 Y</font><br />       A(p3,m,n) = A(a,n)+A(b,m),<font class=\\"jammer\\">4 B) n2 d4 O& r* h- v</font><br />       n--<br /> <span style=\\"display:none\\">; k4 l6 L\\" `' h# v/ c</span>     },<br /> <span style=\\"display:none\\">  m- }3 s* f$ d5 o/ k# i  Y( @+ O</span>     m--<font class=\\"jammer\\">: K. `; r( F' ]/ A3 W- _) j& _\\" w\\" ^</font><br />   },<br /> <span style=\\"display:none\\">: W( ^) @3 n6 z% b9 i6 x</span>   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}<br /> <span style=\\"display:none\\">( ?: e\\" O* Z* d9 O</span> };<br /> <span style=\\"display:none\\">0 l& y\\" S: _6 y) S) n/ n4 g</span>

" c% O( @0 L9 C( U7 X/ o结果：
ndgrid: 3.2001e-002
for1: 1.4999e-002
. C+ ]& `& y9 B0 s, z* @; @) ]for2: 1.86
4 u: E3 x1 t" u0 G& }: Y# Z
0 n6 U" n  W# F/ W0 Y+ N; X1 Z  r+ f

forcal

一段程序的Forcal实现：
9 C; d7 h. j, c8 k1 L% N1 _% `
6 l* e( }3 F- l6 M. {//用C++代码描述为：
: m' F+ V* r! fs=0.0;
for(x=0.0;x<=1.0;x=x+0.0011)
. O6 z( H7 Z: U' l4 {6 Z# f% L{
   for(y=1.0;y<=2.0;y=y+0.0009)
( X4 U5 O" I5 ?, Q$ ~1 C   {
. n1 U8 L4 o5 X+ l" f0 `' p( s     s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   }
}

1、**数组求和函数Sum
3 j) {6 A8 X' y; h5 F) A0 r0 |

2. !using["math","sys"];
   
3.   X+ u6 n0 i9 I1 a7 O2 Q
4. mvar:9 v$ p% ]( J8 G4 W2 i1 A! L. w
   
6. t=clock(),' e( A3 @+ H6 o# h
   
8. oo{7 L, H* e/ Y9 X+ _
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],1 a, Q0 z& ~7 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. * E% g  A+ ?$ r$ I
14. };  ?* N8 ~\\" a' r\\" _; [! k2 c
    
16. [clock()-t]/1000;

!using["math","sys"];<br /> <span style=\\"display:none\\">  X+ u6 n0 i9 I1 a7 O2 Q</span> mvar:<font class=\\"jammer\\">9 v$ p% ]( J8 G4 W2 i1 A! L. w</font><br /> t=clock(),<font class=\\"jammer\\">' e( A3 @+ H6 o# h</font><br /> oo{<font class=\\"jammer\\">7 L, H* e/ Y9 X+ _</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">1 a, Q0 z& ~7 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\\">* E% g  A+ ?$ r$ I</span> };<font class=\\"jammer\\">  ?* N8 ~\\" a' r\\" _; [! k2 c</font><br /> [clock()-t]/1000;

2 v( m$ H4 z0 h- F结果：
1008606.64947441
0.625   //时间
* ?7 }7 }- X/ H
+ a$ V6 e7 R; e2 b& h2 K/ R; r2、求和函数sum
% s) H* u4 L9 c1 w' x) q
; a5 C3 s0 X2 \4 a# m0 t' c

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   4 s9 Y3 [4 E( H# T; b- X' X
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

' x" t7 w' n; X! {结果：
1008606.64947441
6 a" T) n% |  a* ~: r. Q8 C" r0.719   //时间

" j# ^. P1 ?$ y. K- E3、while循环
* c; j3 U; N, I$ z$ C

1. mvar:
   5 k& K6 D8 s$ q6 b+ u6 z
2. t=sys::clock();' K: S, s3 C9 c+ t! t  R
   
3. s=0,x=0, 7 f- \; x! C- V7 l0 {& b8 G
   
4. while{x<=1,  //while循环算法；
   ( `* G9 y\" [: F1 e1 x, W
5.    y=1, 2 W  Q( y  X7 h5 |6 H! }0 x8 }
   
6.    while{y<=2,   s% w  K$ R* R+ ?' H( S
   
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   $ n  B& I. d9 s- l7 @
8.        y=y+0.0009 8 S) ?+ `0 @* k: f
   
9.       },
   & l1 a3 k4 ]+ L; N% x; a
10.    x=x+0.0011
    ! p: b6 e7 [; P5 a5 ?& s$ ~; k/ K
11. },
    , r$ r# `4 y) d: N% y
12. s;  E  a8 ?; N, L& k2 }( S8 I
    
13. [sys::clock()-t]/1000;

复制代码

9 k5 w: ~4 ?& M结果：
1008606.64947441
0.734   //时间
% |% r( E5 ~+ a' X( u6 e

qbist

好深奥！~~~~

forcal

好深奥！~~~~
* X# H/ w9 \- _" u- I6 Lqbist 发表于 2010-10-7 14:56

先了解一下，以备不时之需，有问题可以交流哦，呵呵。

qbist

回复 forcal 的帖子

3 e# _- {# _  n, i3 Q  C
' m+ V2 X5 n* u    嗯！！！

forcal

改进了FcMath中的矩阵乘算法，不知与matlab还有多大差距，朋友们可帮助测一下。

以下是FcMath中的矩阵乘与徐士良算法库XSLSF（普通的C/C++算法）中的矩阵乘的效率比较：
5 R$ ]% K2 L* S% `( U
, A. A( @7 j% q6 I$ |1、FcMath中的矩阵乘
) U' {& `) I9 `, Z1 r- N4 n

2. !using["math","sys"];* w- y2 r: ^9 [. S
   
4. (:a,b,k,t0)=) ^: t9 i4 b% K$ d! f* V- X
   
6. oo{
   
7. 5 n* I5 @- j) ~: X% B5 O& y* d
8.   a=rand[1000,1000], b=rand[1000,1000],; {$ K) a9 T7 M8 t8 |
   
10.   t0=clock(),! s; @6 u. l) k3 _  f5 a1 f. @
    
12.   k=a*b,  //矩阵乘
    
13. $ k& K3 T) S) |! E0 Z: c- V  A
14.   k[1,3:5,9].outm()/ F; V, D- I5 A2 y- j8 I$ Z' M! I
    
16. },
    
17. + X5 Q5 s% c( x
18. [clock()-t0]/1000;
    
19. ! z  |* c  M5 |2 E# w. d% {# p' R0 V( Y

!using["math","sys"];<font class=\\"jammer\\">* w- y2 r: ^9 [. S</font><br /> (:a,b,k,t0)=<font class=\\"jammer\\">) ^: t9 i4 b% K$ d! f* V- X</font><br /> oo{<br /> <span style=\\"display:none\\">5 n* I5 @- j) ~: X% B5 O& y* d</span>   a=rand[1000,1000], b=rand[1000,1000],<font class=\\"jammer\\">; {$ K) a9 T7 M8 t8 |</font><br />   t0=clock(),<font class=\\"jammer\\">! s; @6 u. l) k3 _  f5 a1 f. @</font><br />   k=a*b,  //矩阵乘<br /> <span style=\\"display:none\\">$ k& K3 T) S) |! E0 Z: c- V  A</span>   k[1,3:5,9].outm()<font class=\\"jammer\\">/ F; V, D- I5 A2 y- j8 I$ Z' M! I</font><br /> },<br /> <span style=\\"display:none\\">+ X5 Q5 s% c( x</span> [clock()-t0]/1000;<br /> <span style=\\"display:none\\">! z  |* c  M5 |2 E# w. d% {# p' R0 V( Y</span>

结果：
0 Z0 n' E8 D2 W1 u9 x. _( K+ V

1.         238.447        247.837        247.065        248.105        247.0587 k. u2 M7 t\" O1 k+ E# ]$ o
   
2.         244.123        249.925        247.553        243.981        250.016& h# G& Y( A) @0 A& R; O\" m4 u0 A
   
3.         236.387        252.025        245.651        248.866        248.866+ e/ i( f1 ?5 l8 _
   
4. 2.219 秒
     @4 `$ z% b# z  Q

复制代码

2、XSLSF（普通的C/C++算法）中的矩阵乘
% U& j& f# f! h8 C

2. !using["math","sys","XSLSF"];' e; l: K# o1 e# _2 P2 R
   
4. (:a,b,k,t0)=
   
5. * c2 A& X* N4 I- P) G
6. oo{
   
7. % t  p  K# ~% g1 c7 A4 r
8.   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],) e\\" W$ q+ _1 S  ?8 _' v+ @7 J. z5 J
   
10.   t0=clock(),
    
11. . Q+ Q% n) \  p5 Y
12.   rmul[k:a,b],  //矩阵乘8 f( J# d. x# m) Z+ Z; d2 r
    
14.   k[1,3:5,9].outm()
    
15. 8 B9 J7 F7 }3 k% ]7 G( o& j- ^
16. },
    
17. , P! T0 Y; b$ N
18. [clock()-t0]/1000;2 M$ C3 R& z- }* J) |' e
    

!using["math","sys","XSLSF"];<font class=\\"jammer\\">' e; l: K# o1 e# _2 P2 R</font><br /> (:a,b,k,t0)=<br /> <span style=\\"display:none\\">* c2 A& X* N4 I- P) G</span> oo{<br /> <span style=\\"display:none\\">% t  p  K# ~% g1 c7 A4 r</span>   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],<font class=\\"jammer\\">) e\\" W$ q+ _1 S  ?8 _' v+ @7 J. z5 J</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">. Q+ Q% n) \  p5 Y</span>   rmul[k:a,b],  //矩阵乘<font class=\\"jammer\\">8 f( J# d. x# m) Z+ Z; d2 r</font><br />   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">8 B9 J7 F7 }3 k% ]7 G( o& j- ^</span> },<br /> <span style=\\"display:none\\">, P! T0 Y; b$ N</span> [clock()-t0]/1000;<font class=\\"jammer\\">2 M$ C3 R& z- }* J) |' e</font><br />

结果：

1.         262.121        247.583        260.529        259.548        258.328, r, ~6 N; b( n1 [; ]# e
   
2.         255.413        246.563        254.356        250.548        251.5094 u- a! C* l5 e$ N8 T! V& j
   
3.         256.152        247.725        259.444        250.827        249.816
   ; r* G8 v6 e' Z- H8 x* `( i* C/ g
4. 10.563 秒
   0 G# L! K9 L& f\" r

复制代码

8 v' @! n# B& C1 p