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的数学爱好者的极佳的练手工具，同时也期望对一般的数值计算用户提供越来越多的方便。

: X$ T; f& \* F

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

9 |2 ]7 r: p' R. G; Z) ~$ o1 e6 k0 W

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

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

forcal

例子1代码：

2. !using["math"];
   
3. % Y0 {0 k3 S& X9 L\\" v3 ^6 _
4. mvar:# p2 n3 ?: e# H* A1 c+ B
   
6. oo{                      //一般在oo函数中调用FcMath函数
   
7. 1 B4 C0 P  C% J: v# w
8.   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化
   
9. # {5 [! F1 q/ i8 L7 b; A7 S0 q
10.   a.outm(),              //输出矩阵a/ l* ^6 p/ q; I6 K- p4 a* U$ e
    
12.   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出; ?\\" V3 q\\" |( u% H+ Y( O* m
    
14.   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出
    
15. 4 N  S) j4 s) V
16.   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出
    
17. \\" G7 Q+ D, G/ o4 ]$ k1 W1 W
18. };# L! u/ P1 L0 O+ J' y6 [
    

!using["math"];<br /> <span style=\\"display:none\\">% Y0 {0 k3 S& X9 L\\" v3 ^6 _</span> mvar:<font class=\\"jammer\\"># p2 n3 ?: e# H* A1 c+ B</font><br /> oo{                      //一般在oo函数中调用FcMath函数<br /> <span style=\\"display:none\\">1 B4 C0 P  C% J: v# w</span>   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化<br /> <span style=\\"display:none\\"># {5 [! F1 q/ i8 L7 b; A7 S0 q</span>   a.outm(),              //输出矩阵a<font class=\\"jammer\\">/ l* ^6 p/ q; I6 K- p4 a* U$ e</font><br />   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">; ?\\" V3 q\\" |( u% H+ Y( O* m</font><br />   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出<br /> <span style=\\"display:none\\">4 N  S) j4 s) V</span>   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出<br /> <span style=\\"display:none\\">\\" G7 Q+ D, G/ o4 ]$ k1 W1 W</span> };<font class=\\"jammer\\"># L! u/ P1 L0 O+ J' y6 [</font><br />

结果：

1.        0.211319   4.91638e-002       0.144638       0.153259       0.852615- S* L' h0 u8 x\" U& h  v8 h# z+ ?! e
   
2.        0.630646       0.927048       0.440308       0.162857       0.556854: W  Y& E  [2 i2 j\" b\" ^! C1 N, c
   
3.         0.43309        0.34552       0.563919       0.937164       0.209641
   ) h' ], B; _! {# o- y' ~* u
4.        0.603271       0.727676       0.130951   5.35736e-002       0.197937
   3 O2 {7 n4 X8 X8 H. F. H! |+ }
5.        0.576004       0.747589   1.17645e-002       0.363892       0.280777
   # ?' k) }* o* ^4 }- V) X
6.        0.646454       0.381088        0.58551        0.26387        0.93692
   # P+ T; a: Q8 K; e( s
7. 
   & t1 T( J6 x6 ?! v! f* a9 x
8.        0.153259
   ' Q% Z' J\" q$ I
9.        0.162857. A3 @. R5 _3 k/ Y
   
10.        0.937164
    3 i6 T) E( x& S' I* G6 l  s! |9 p
11.    5.35736e-002  l2 H; X7 [' K+ K4 H
    
12.        0.363892
    ! s6 [\" F2 X+ t
13.         0.26387  a( r* k8 F& F\" b7 x3 V7 S5 {( x( w
    
14. 
    % A\" S3 ]/ W8 M( s, s\" O
15.        0.603271       0.727676       0.130951   5.35736e-002       0.197937! }2 z' l% c7 `! O9 ^' U
    
16. 
    9 q$ ~6 W\" L4 a
17.        0.130951   5.35736e-002
    # @3 F6 Q5 y0 h. w
18.    1.17645e-002       0.363892
    5 n# q3 \2 b, W\" d2 ^& v+ t
19.         0.58551        0.26387
    1 ?* G\" u' D  t) I& f: e
20. 0 _! a( I& Y4 Q
    

复制代码

例子2代码：
* b1 f# d7 r% Y6 e& M+ w- M: G
$ |. X9 H/ d. Q2 j5 D# F( u! Y

2. f(x1,x2,x3,y1,y2,y3)=      //函数定义( y; |1 j, G$ y3 g8 K
   
4. {
   
5. % n: u, b/ d9 _! Q: y1 |8 ~! {
6.     y1=x1*x1+x2*x2+x3*x3-1.0,7 j6 u) Y# M8 F& s5 |
   
8.     y2=2.0*x1*x1+x2*x2-4.0*x3,
   
9. 2 y4 b8 \; C; }. u+ x
10.     y3=3.0*x1*x1-4.0*x2+x3*x3
    
11. , v* d% r. j. S& Y+ M0 x
12. };
    
13. , `, m5 I1 Z) C3 t& y
14. !using["math","sys"];9 H7 T\\" N) ^, e5 g, J
    
16. mvar:/ ]  u, @+ b' V' O) U( u; x
    
18. oo{3 m) h5 E# C5 L6 x
    
20.   x=array(3),
    
21. ) b; ]  L5 [+ s# S( p& b
22.   x.SA[0 : 1,1,1],       //设置初值为1,1,1
    
23. & K0 L* A/ x0 y  J
24.   i=netn[HFor("f"),x],   //拟牛顿法解方程
    
25. 0 ]! ]7 l' D3 j5 {1 {+ T$ S
26.   x.outm(),              //输出结果
    
27. * B+ _/ P) l# W8 S
28.   i                      //返回迭代次数: l3 G4 X7 A8 h8 o/ X9 w
    
30. };$ N8 N1 q% F- [: C- a0 `; V# l5 W
    

f(x1,x2,x3,y1,y2,y3)=      //函数定义<font class=\\"jammer\\">( y; |1 j, G$ y3 g8 K</font><br /> {<br /> <span style=\\"display:none\\">% n: u, b/ d9 _! Q: y1 |8 ~! {</span>     y1=x1*x1+x2*x2+x3*x3-1.0,<font class=\\"jammer\\">7 j6 u) Y# M8 F& s5 |</font><br />     y2=2.0*x1*x1+x2*x2-4.0*x3,<br /> <span style=\\"display:none\\">2 y4 b8 \; C; }. u+ x</span>     y3=3.0*x1*x1-4.0*x2+x3*x3<br /> <span style=\\"display:none\\">, v* d% r. j. S& Y+ M0 x</span> };<br /> <span style=\\"display:none\\">, `, m5 I1 Z) C3 t& y</span> !using["math","sys"];<font class=\\"jammer\\">9 H7 T\\" N) ^, e5 g, J</font><br /> mvar:<font class=\\"jammer\\">/ ]  u, @+ b' V' O) U( u; x</font><br /> oo{<font class=\\"jammer\\">3 m) h5 E# C5 L6 x</font><br />   x=array(3),<br /> <span style=\\"display:none\\">) b; ]  L5 [+ s# S( p& b</span>   x.SA[0 : 1,1,1],       //设置初值为1,1,1<br /> <span style=\\"display:none\\">& K0 L* A/ x0 y  J</span>   i=netn[HFor("f"),x],   //拟牛顿法解方程<br /> <span style=\\"display:none\\">0 ]! ]7 l' D3 j5 {1 {+ T$ S</span>   x.outm(),              //输出结果<br /> <span style=\\"display:none\\">* B+ _/ P) l# W8 S</span>   i                      //返回迭代次数<font class=\\"jammer\\">: l3 G4 X7 A8 h8 o/ X9 w</font><br /> };<font class=\\"jammer\\">$ N8 N1 q% F- [: C- a0 `; V# l5 W</font><br />

结果：
  0.785197       0.496611       0.369923
" W( r$ B0 G+ k# M; y

forcal

效率测试：
/ A3 V6 ~. q+ k- Y1 h1 Zsimwe的网友lin2009 的matlab代码：

1. clear all4 ~6 F( d8 X  S9 [2 ]! U+ y
   
2. clc7 C- K8 o# D9 Y9 a
   
3. tic' `2 U8 N/ y8 y6 b9 J& h$ Z* O
   
4. k = zeros(5,5); % //生成5×5全0矩阵6 S: h7 r7 @' M$ |
   
5. % 循环计算以下程序段100000次：% t2 C3 f9 `! h# T4 C
   
6. for m = 1:1000007 R2 Q8 H5 k. r3 T
   
7.     a = rand(5,7);  s+ W% G/ D0 U9 C) B
   
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化) W2 h+ O, l; H, N0 e7 C. U  `9 K
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);$ @* p/ n9 l5 w1 G6 R! j
   
10. end
    0 C. P\" L& R4 Y: e+ W% ~
11. k
    9 K9 m, W2 b! ~, E* X
12. toc
    % z; H. A7 g3 B) U; u- [0 S; D

复制代码

Forcal代码：
# i( y7 \# u! b' J, G' Y5 G3 r7 X
运行稍快的代码，比matlab约快10％吧？

2. !using["math","sys"];6 |5 t% [6 p% |5 s+ O
   
4. mvar:
   
5. ! I  t- {9 i# m: L1 `
6. t0=clock(),
   
7. & O9 z0 X9 p, q6 I
8. oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为01 P  |6 _, j7 R# U6 t
   
10. i=0,(i<1000 00).while{ //循环计算1000 00次1 U  Q( t# V( p8 d
    
12.   oo{
    
13. 7 g8 E* z2 T3 R* g
14.     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
    
15. ( _! R8 L. l+ w3 S- w5 Q  ]' c1 g
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. 7 L: }+ }- d( K: n6 v- k
18.   },
    
19. \\" m\\" `: K* b3 R
20.   i++
    
21. 4 V; |2 m- P1 m; a6 X
22. },
    
23. 2 J: Y& y# k# C; x
24. k.outm(),             //输出矩阵k，然后销毁k
    
25. & O' n0 D9 n  a2 I
26. [clock()-t0]/1000;    //得到计算时间，秒

!using["math","sys"];<font class=\\"jammer\\">6 |5 t% [6 p% |5 s+ O</font><br /> mvar:<br /> <span style=\\"display:none\\">! I  t- {9 i# m: L1 `</span> t0=clock(),<br /> <span style=\\"display:none\\">& O9 z0 X9 p, q6 I</span> oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0<font class=\\"jammer\\">1 P  |6 _, j7 R# U6 t</font><br /> i=0,(i<1000 00).while{ //循环计算1000 00次<font class=\\"jammer\\">1 U  Q( t# V( p8 d</font><br />   oo{<br /> <span style=\\"display:none\\">7 g8 E* z2 T3 R* g</span>     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化<br /> <span style=\\"display:none\\">( _! R8 L. l+ w3 S- w5 Q  ]' c1 g</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\\">7 L: }+ }- d( K: n6 v- k</span>   },<br /> <span style=\\"display:none\\">\\" m\\" `: K* b3 R</span>   i++<br /> <span style=\\"display:none\\">4 V; |2 m- P1 m; a6 X</span> },<br /> <span style=\\"display:none\\">2 J: Y& y# k# C; x</span> k.outm(),             //输出矩阵k，然后销毁k<br /> <span style=\\"display:none\\">& O' n0 D9 n  a2 I</span> [clock()-t0]/1000;    //得到计算时间，秒

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

% O; s+ a, O0 p# O9 c1 d$ s8 R比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];
   
3. . a. N1 y4 \1 ~/ u5 P
4. (:t0,k,i,a,b)=
   
5. ! n- @& ?, w/ D& D7 [! s3 g  x0 X
6. {. W: Y) E4 c) A0 F& u4 X& T
   
8.   t0=clock(),7 E! a6 Z0 b$ i  z\\" j# s
   
10.   k=zeros[5,5],
    
11. / B, |. M' U7 ?* I
12.   i=0,(i<1000 00).while{
    
13. \\" K# U- n2 p- L0 T2 a9 f' P: S\\" |
14.     oo{
    
15. 2 i8 y6 s; o. E0 E, |. `* M
16.       a=rand[5,7], b=rand[7,5],
    
17. # P6 P\\" ~5 o) x3 ~
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg), P1 y: ^4 v7 ?+ g$ u8 q) S
    
20.     },7 y$ ?: ?8 L8 ~; i
    
22.     i++
    
23. , B8 {# ^# h- m: n
24.   },7 ?7 M$ @+ I/ b
    
26.   k.outm().delete(),
    
27. ! h. o' N# W2 `
28.   [clock()-t0]/1000
    
29.   b0 m8 A9 o/ \( k- a+ m9 c
30. };

!using["math","sys"];<br /> <span style=\\"display:none\\">. a. N1 y4 \1 ~/ u5 P</span> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">! n- @& ?, w/ D& D7 [! s3 g  x0 X</span> {<font class=\\"jammer\\">. W: Y) E4 c) A0 F& u4 X& T</font><br />   t0=clock(),<font class=\\"jammer\\">7 E! a6 Z0 b$ i  z\\" j# s</font><br />   k=zeros[5,5],<br /> <span style=\\"display:none\\">/ B, |. M' U7 ?* I</span>   i=0,(i<1000 00).while{<br /> <span style=\\"display:none\\">\\" K# U- n2 p- L0 T2 a9 f' P: S\\" |</span>     oo{<br /> <span style=\\"display:none\\">2 i8 y6 s; o. E0 E, |. `* M</span>       a=rand[5,7], b=rand[7,5],<br /> <span style=\\"display:none\\"># P6 P\\" ~5 o) x3 ~</span>       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<font class=\\"jammer\\">, P1 y: ^4 v7 ?+ g$ u8 q) S</font><br />     },<font class=\\"jammer\\">7 y$ ?: ?8 L8 ~; i</font><br />     i++<br /> <span style=\\"display:none\\">, B8 {# ^# h- m: n</span>   },<font class=\\"jammer\\">7 ?7 M$ @+ I/ b</font><br />   k.outm().delete(),<br /> <span style=\\"display:none\\">! h. o' N# W2 `</span>   [clock()-t0]/1000<br /> <span style=\\"display:none\\">  b0 m8 A9 o/ \( k- a+ m9 c</span> };

$ c: R% p; v# y/ F- |. ^& A$ W在我的电脑上运行时间为3.579秒。
4 [6 w5 Y* ]4 P5 L+ O3 ^
该例子的理论结果是每个元素均为275000。
: v  Q( a* }/ P2 _1 {
我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。
( v6 M' W9 G1 t% {4 a- d! \

forcal

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

2. !using["math"];% R7 s$ c+ k4 x7 r\\" T  [% A
   
4. mvar:( b+ e9 k- A! t+ ?3 p
   
6. oo{
   
7. ' w0 P; ^' A' `1 }6 @
8.   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
   
9. ' u# v9 |) e0 ?& V& Q
10.   a1.outm[5,1,1],! ~3 Y/ W* z7 a3 |# N4 c5 P
    
12.   a2.outm[5,1,1],
    
13. 9 C/ g8 |5 F0 H1 e! I  b
14.   a3.outm[5,1,1],  D! _# h2 y8 ^5 I- ^3 S# G# j
    
16.   a4.outm[5,1,1],
    
17. 3 O. Q7 L( b. l2 A' f) Y
18.   a=a1+a2+a3+a4,
    
19. ' U% }6 ?& p) ]& _/ n6 U
20.   a.outm[5,1,1],
    
21. 7 v+ e9 j/ b4 |; r. A; b
22.   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]$ P\\" o! F( W8 l+ ]. u( q8 F: ^
    
24. };6 @8 p3 F  Z- U! r- P
    

!using["math"];<font class=\\"jammer\\">% R7 s$ c+ k4 x7 r\\" T  [% A</font><br /> mvar:<font class=\\"jammer\\">( b+ e9 k- A! t+ ?3 p</font><br /> oo{<br /> <span style=\\"display:none\\">' w0 P; ^' A' `1 }6 @</span>   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],<br /> <span style=\\"display:none\\">' u# v9 |) e0 ?& V& Q</span>   a1.outm[5,1,1],<font class=\\"jammer\\">! ~3 Y/ W* z7 a3 |# N4 c5 P</font><br />   a2.outm[5,1,1],<br /> <span style=\\"display:none\\">9 C/ g8 |5 F0 H1 e! I  b</span>   a3.outm[5,1,1],<font class=\\"jammer\\">  D! _# h2 y8 ^5 I- ^3 S# G# j</font><br />   a4.outm[5,1,1],<br /> <span style=\\"display:none\\">3 O. Q7 L( b. l2 A' f) Y</span>   a=a1+a2+a3+a4,<br /> <span style=\\"display:none\\">' U% }6 ?& p) ]& _/ n6 U</span>   a.outm[5,1,1],<br /> <span style=\\"display:none\\">7 v+ e9 j/ b4 |; r. A; b</span>   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]<font class=\\"jammer\\">$ P\\" o! F( W8 l+ ]. u( q8 F: ^</font><br /> };<font class=\\"jammer\\">6 @8 p3 F  Z- U! r- P</font><br />

3 B! X+ H3 K* q4 B: l说明：
6 U4 G8 d) g- O+ X: Y" P' j+ rlinspace(8,12,5)：生成一维数组，共5个元素8～12
a1.outm[5,1,1]：输出**数组a1，连下标一起输出
. C7 D4 t6 E, ]8 ySum[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)。
4 V4 e  f. a. h- C" g/ Z
结果（最终求和结果是1320）：

(0,0,0,*)              1.            1.            1.            1.            1.
(0,0,1,*)              1.            1.            1.            1.            1.
# C" Q7 A/ W1 H' W; \(0,1,0,*)              1.            1.            1.            1.            1.
(0,1,1,*)              1.            1.            1.            1.            1.
(0,2,0,*)              1.            1.            1.            1.            1.
8 [  B" ^$ g3 y" B& U% [(0,2,1,*)              1.            1.            1.            1.            1.
(1,0,0,*)              2.            2.            2.            2.            2.
( V7 L. U1 N/ K. e$ r! p(1,0,1,*)              2.            2.            2.            2.            2.
(1,1,0,*)              2.            2.            2.            2.            2.
(1,1,1,*)              2.            2.            2.            2.            2.
- y$ d1 p; L: i  Z+ h(1,2,0,*)              2.            2.            2.            2.            2.
(1,2,1,*)              2.            2.            2.            2.            2.
1 G4 M9 ]' I4 N/ `
4 p! |  z3 G8 b/ E. A3 ~, q(0,0,0,*)              3.            3.            3.            3.            3.
(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.
(0,2,1,*)              5.            5.            5.            5.            5.
(1,0,0,*)              3.            3.            3.            3.            3.
(1,0,1,*)              3.            3.            3.            3.            3.
, Y+ f; [# R# f% |+ t3 s4 N, U(1,1,0,*)              4.            4.            4.            4.            4.
: v" F) a7 x3 Z(1,1,1,*)              4.            4.            4.            4.            4.
(1,2,0,*)              5.            5.            5.            5.            5.
3 `5 X& ]3 l& i) a2 Q8 {: t% x(1,2,1,*)              5.            5.            5.            5.            5.

6 T8 a7 r6 W: U* K+ B! e3 O(0,0,0,*)              6.            6.            6.            6.            6.
8 Q0 G8 I  ?8 C. V: b- j% L(0,0,1,*)              7.            7.            7.            7.            7.
(0,1,0,*)              6.            6.            6.            6.            6.
(0,1,1,*)              7.            7.            7.            7.            7.
* ?5 q* C/ o2 c(0,2,0,*)              6.            6.            6.            6.            6.
! p) l/ M* X5 u" k(0,2,1,*)              7.            7.            7.            7.            7.
- ?. N  u5 _" I- E( }(1,0,0,*)              6.            6.            6.            6.            6.
(1,0,1,*)              7.            7.            7.            7.            7.
# D# U  M/ i. X/ t4 w2 t(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.
* g" w/ C: i+ K5 \
$ \7 N' T: V+ l8 R/ K(0,0,0,*)              8.            9.           10.           11.           12.
# M- j9 ~( }1 }9 x" U3 \(0,0,1,*)              8.            9.           10.           11.           12.
(0,1,0,*)              8.            9.           10.           11.           12.
. S5 l; B% G3 \8 r7 ^( \(0,1,1,*)              8.            9.           10.           11.           12.
(0,2,0,*)              8.            9.           10.           11.           12.
(0,2,1,*)              8.            9.           10.           11.           12.
& n# e; L3 g0 ^3 f, ]+ S+ k4 ~(1,0,0,*)              8.            9.           10.           11.           12.
. B. B7 H0 o( [(1,0,1,*)              8.            9.           10.           11.           12.
(1,1,0,*)              8.            9.           10.           11.           12.
& e, `' Y9 y7 P(1,1,1,*)              8.            9.           10.           11.           12.
(1,2,0,*)              8.            9.           10.           11.           12.
(1,2,1,*)              8.            9.           10.           11.           12.

(0,0,0,*)             18.           19.           20.           21.           22.
(0,0,1,*)             19.           20.           21.           22.           23.
5 X" t: S9 C" w4 B(0,1,0,*)             19.           20.           21.           22.           23.
: b9 }6 @8 \0 T/ g+ ~( O: ?(0,1,1,*)             20.           21.           22.           23.           24.
(0,2,0,*)             20.           21.           22.           23.           24.
) h9 W6 _7 g$ c; A) k(0,2,1,*)             21.           22.           23.           24.           25.
(1,0,0,*)             19.           20.           21.           22.           23.
(1,0,1,*)             20.           21.           22.           23.           24.
(1,1,0,*)             20.           21.           22.           23.           24.
(1,1,1,*)             21.           22.           23.           24.           25.
9 F+ @, P* W! b(1,2,0,*)             21.           22.           23.           24.           25.
- g- |4 E4 p: t( B3 ~(1,2,1,*)             22.           23.           24.           25.           26.

(0,0,*)            100.          105.
1 I6 D: \' ?2 |' O(0,1,*)            105.          110.
& D" u2 C; g; x3 w- A(0,2,*)            110.          115.
2 G) C# ?! O+ m4 _2 ]/ B" E+ c(1,0,*)            105.          110.
9 U0 M+ s9 }9 l  I6 \(1,1,*)            110.          115.
(1,2,*)            115.          120.
6 b$ j# s+ e+ E4 }  k2 I
7 O4 J- U/ g+ }, U4 ]# D  i(0,*)            205.          215.          225.
& x' E" l; r( j' y3 E(1,*)            215.          225.          235.
# a) t/ ^( n7 V/ l9 {8 |; h$ \
5 f2 Q4 L/ m# b. R: M; u(0,*)            645.
(1,*)            675.
: f6 b9 ]' n8 R+ B$ Z& k
4 Q5 K8 y" e# _1320.
  b( K+ B+ J5 s  d* _/ K6 ]2 ?$ ?1 c

forcal

在matlab中，纯for循环速度最慢。而一半for循环+一半向量化的速度最快，Forcal中也是如此：

2. !using["math","sys"];8 s6 D/ X+ l4 ]8 P* m; ~
   
4. mvar:$ k: \$ D  ^( c, E* {
   
6. (:p1,p2,p3,a,b)=
   
7. 7 S6 P: C, z2 j, ]  H+ g
8. {
   
9.   W# @' |% A& \# u
10.   oo{
    
11. ' ]8 X& j3 V( N8 R( W
12.     a=array[1000].rand(),
    
13. / |, I4 }' R# z# y' O
14.     b=array[1000].rand(),( A' g+ ^, ?/ g# t7 }3 {
    
16.     p1=array[1000,1000],6 J, ~4 A0 A# X; r) `
    
18.     p2=array[1000,1000],
    
19. . d% u7 t0 ], b& O. [  w! D
20.     p3=array[1000,1000],. i* A5 x6 F4 o! Z6 O7 y' C; u
    
22.     t0=clock(),
    
23. * c6 l- O, p& a& O% e# C% h
24.     ndgrid(a,b,&A,&B),
    
25. $ ?* V% v9 O. F8 N( ?4 E
26.     p1.=A+B; u0 ?: F7 o! s  V8 _! @. Q
    
28.   },
    
29. & O0 j  f1 J  ?) q, a9 j
30.   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},9 G9 @  a% [( u! v
    
32.   lena=FCDLen(a),/ m1 f/ K5 M7 s  q0 X
    
34.   lenb=FCDLen(b),$ M5 e- `! e$ N: V
    
36.   t0=clock(),( U0 A0 ~2 }  j
    
38.   m = lenb-1, (m>=0).while{
    
39. $ Q\\" b$ a$ d! e4 b% i; Z2 s! P
40.     oo{p2(m,neg) = a+rn[b(m)]},
    
41. 7 d3 V$ `' c; L; C1 [. Q! C' w# c& v
42.     m--
    
43. 4 F4 \9 B1 R$ c  ]& M  U0 D. ]
44.   },
    
45. , l\\" _* r0 u. Z4 [' K+ H# {. K
46.   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},
    
47. 8 P! S. C2 e. v- L% @
48.   t0=clock(),6 B/ c6 e. i, B+ P- M
    
50.   m = lenb-1, (m>=0).while{
    
51. - z\\" p$ S8 p, B; V; W9 X
52.     n = lena-1, (n>=0).while{\\" \3 e: ?/ |. w2 w/ l' u: r) I
    
54.       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些# r- i: w0 R, R\\" c- k+ z
    
56.       A(p3,m,n) = A(a,n)+A(b,m),
    
57. ) }  _9 J0 g- C+ d\\" u# Y* ^
58.       n--
    
59. 1 y- R6 a/ P3 V! j3 s
60.     },# m; n\\" Y1 @  ?- W  ?
    
62.     m--+ ]2 n0 f! c- h& V9 x9 l5 P\\" a5 c
    
64.   },3 H$ V\\" {9 z- e) I\\" d/ o
    
66.   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
    
67. 6 H) N' C. Z4 S' \: P9 `- A# e; n
68. };/ v9 f2 g( d, a4 N+ i
    

!using["math","sys"];<font class=\\"jammer\\">8 s6 D/ X+ l4 ]8 P* m; ~</font><br /> mvar:<font class=\\"jammer\\">$ k: \$ D  ^( c, E* {</font><br /> (:p1,p2,p3,a,b)=<br /> <span style=\\"display:none\\">7 S6 P: C, z2 j, ]  H+ g</span> {<br /> <span style=\\"display:none\\">  W# @' |% A& \# u</span>   oo{<br /> <span style=\\"display:none\\">' ]8 X& j3 V( N8 R( W</span>     a=array[1000].rand(),<br /> <span style=\\"display:none\\">/ |, I4 }' R# z# y' O</span>     b=array[1000].rand(),<font class=\\"jammer\\">( A' g+ ^, ?/ g# t7 }3 {</font><br />     p1=array[1000,1000],<font class=\\"jammer\\">6 J, ~4 A0 A# X; r) `</font><br />     p2=array[1000,1000],<br /> <span style=\\"display:none\\">. d% u7 t0 ], b& O. [  w! D</span>     p3=array[1000,1000],<font class=\\"jammer\\">. i* A5 x6 F4 o! Z6 O7 y' C; u</font><br />     t0=clock(),<br /> <span style=\\"display:none\\">* c6 l- O, p& a& O% e# C% h</span>     ndgrid(a,b,&A,&B),<br /> <span style=\\"display:none\\">$ ?* V% v9 O. F8 N( ?4 E</span>     p1.=A+B<font class=\\"jammer\\">; u0 ?: F7 o! s  V8 _! @. Q</font><br />   },<br /> <span style=\\"display:none\\">& O0 j  f1 J  ?) q, a9 j</span>   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},<font class=\\"jammer\\">9 G9 @  a% [( u! v</font><br />   lena=FCDLen(a),<font class=\\"jammer\\">/ m1 f/ K5 M7 s  q0 X</font><br />   lenb=FCDLen(b),<font class=\\"jammer\\">$ M5 e- `! e$ N: V</font><br />   t0=clock(),<font class=\\"jammer\\">( U0 A0 ~2 }  j</font><br />   m = lenb-1, (m>=0).while{<br /> <span style=\\"display:none\\">$ Q\\" b$ a$ d! e4 b% i; Z2 s! P</span>     oo{p2(m,neg) = a+rn[b(m)]},<br /> <span style=\\"display:none\\">7 d3 V$ `' c; L; C1 [. Q! C' w# c& v</span>     m--<br /> <span style=\\"display:none\\">4 F4 \9 B1 R$ c  ]& M  U0 D. ]</span>   },<br /> <span style=\\"display:none\\">, l\\" _* r0 u. Z4 [' K+ H# {. K</span>   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},<br /> <span style=\\"display:none\\">8 P! S. C2 e. v- L% @</span>   t0=clock(),<font class=\\"jammer\\">6 B/ c6 e. i, B+ P- M</font><br />   m = lenb-1, (m>=0).while{<br /> <span style=\\"display:none\\">- z\\" p$ S8 p, B; V; W9 X</span>     n = lena-1, (n>=0).while{<font class=\\"jammer\\">\\" \3 e: ?/ |. w2 w/ l' u: r) I</font><br />       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些<font class=\\"jammer\\"># r- i: w0 R, R\\" c- k+ z</font><br />       A(p3,m,n) = A(a,n)+A(b,m),<br /> <span style=\\"display:none\\">) }  _9 J0 g- C+ d\\" u# Y* ^</span>       n--<br /> <span style=\\"display:none\\">1 y- R6 a/ P3 V! j3 s</span>     },<font class=\\"jammer\\"># m; n\\" Y1 @  ?- W  ?</font><br />     m--<font class=\\"jammer\\">+ ]2 n0 f! c- h& V9 x9 l5 P\\" a5 c</font><br />   },<font class=\\"jammer\\">3 H$ V\\" {9 z- e) I\\" d/ o</font><br />   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}<br /> <span style=\\"display:none\\">6 H) N' C. Z4 S' \: P9 `- A# e; n</span> };<font class=\\"jammer\\">/ v9 f2 g( d, a4 N+ i</font><br />

结果：
ndgrid: 3.2001e-002
for1: 1.4999e-002
for2: 1.86
* W7 t6 N' l2 |) j4 ]

forcal

一段程序的Forcal实现：
: @. {* ^5 u0 u
//用C++代码描述为：
4 p( T9 f3 ~- c; }3 `5 B- G( gs=0.0;
4 ]' F; y0 \+ H! R6 `" Y3 E3 Yfor(x=0.0;x<=1.0;x=x+0.0011)
- ~4 z% U' H) p' x- S* {{
+ G# l3 \6 D2 n: g/ ?   for(y=1.0;y<=2.0;y=y+0.0009)
   {
     s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
+ V: j, }# }% c; o: I' i9 E0 Z+ P   }
}

4 f6 W; X: J* i# [# ?  |6 t% U1、**数组求和函数Sum
9 Z0 b+ N$ ]- x3 P( W5 ]. V
! A: |2 x- y; \3 h$ O

2. !using["math","sys"];' U+ ~4 G' }/ U% @3 K
   
4. mvar:7 ~/ K$ w; Q) Z9 K7 G
   
6. t=clock(),8 {, X% @\\" J& R2 C; W
   
8. oo{
   
9. & C, O. E. `$ W! ~1 Q
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],% P0 ~3 n7 c8 N
    
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]: Q3 h) ?1 [2 l3 V5 P/ }+ N0 v# M\\" N2 C
    
14. };  Y, J# `# G3 P3 z; M  T
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">' U+ ~4 G' }/ U% @3 K</font><br /> mvar:<font class=\\"jammer\\">7 ~/ K$ w; Q) Z9 K7 G</font><br /> t=clock(),<font class=\\"jammer\\">8 {, X% @\\" J& R2 C; W</font><br /> oo{<br /> <span style=\\"display:none\\">& C, O. E. `$ W! ~1 Q</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">% P0 ~3 n7 c8 N</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\\">: Q3 h) ?1 [2 l3 V5 P/ }+ N0 v# M\\" N2 C</font><br /> };<font class=\\"jammer\\">  Y, J# `# G3 P3 z; M  T</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
0.625   //时间
1 r$ H  L7 |' `# m2 G' y
  K3 D5 z$ h" }. `: H, X5 h2、求和函数sum
+ j7 R: U6 }2 H+ ]2 z

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   & {( X# u2 C3 u0 b
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
0 V- L6 ^0 a2 f& l# v1008606.64947441
0.719   //时间
' {4 O9 }, k9 Y* p- i5 M  Z
( ]6 D& [. W% \4 H1 k3 e3、while循环
/ q8 N0 ?% w: k7 a2 q5 [# s* E
" m! Z7 @$ I% W: n. p4 Q% I1 ]6 S

1. mvar:
   9 O: L! H8 j6 ~+ w) w6 p, U
2. t=sys::clock();4 p+ w9 {! O4 n
   
3. s=0,x=0, # J6 ]9 x9 c) J& \
   
4. while{x<=1,  //while循环算法；
   + u$ @9 j/ H8 n$ R% r& |, Y
5.    y=1,
   0 i' K( I* K! o$ _. w4 ~& U
6.    while{y<=2,
   6 {; T+ @! x: W  e& `
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   ; p; i: l4 C5 n% w8 i
8.        y=y+0.0009
   3 S& e) K5 ]\" e/ L2 {
9.       }, 3 S7 }$ [7 T, ^7 P; N# B9 j( M
   
10.    x=x+0.0011 / }6 x; e+ s& o) K; `/ t0 N/ Z# A
    
11. },
    ( ?8 Q% g, y) |7 v2 K6 V. V
12. s;6 u4 p, G% Y\" x+ ^
    
13. [sys::clock()-t]/1000;

复制代码

3 W1 v3 _$ y' F8 P结果：
" G! i, ]/ q- `1008606.64947441
0.734   //时间
( H% |3 M1 ^$ R2 B, c5 X7 B

qbist

好深奥！~~~~

forcal

好深奥！~~~~
, z  f5 f, Y1 Z( }# @; ]qbist 发表于 2010-10-7 14:56

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

qbist

回复 forcal 的帖子


    嗯！！！

forcal

改进了FcMath中的矩阵乘算法，不知与matlab还有多大差距，朋友们可帮助测一下。
- |  d$ i3 c+ c3 O/ w
% w8 ^: G+ B- X: b以下是FcMath中的矩阵乘与徐士良算法库XSLSF（普通的C/C++算法）中的矩阵乘的效率比较：

) _1 r% X2 h' w" r% w4 K: f1、FcMath中的矩阵乘

2. !using["math","sys"];0 a6 k/ k\\" V0 z. k! k
   
4. (:a,b,k,t0)=
   
5. 2 s! P; V8 [* w8 s
6. oo{3 V9 O\\" q7 I. Y2 N6 Q- ?, J& j
   
8.   a=rand[1000,1000], b=rand[1000,1000],
   
9. & t4 i: ^4 g: W6 R
10.   t0=clock(),
    
11. & X4 G% p. k+ t! T
12.   k=a*b,  //矩阵乘% n+ O9 I$ }4 y+ Z
    
14.   k[1,3:5,9].outm()
    
15. 0 X8 C# }+ d' J
16. },7 M# {4 x; D6 C
    
18. [clock()-t0]/1000;& l% Q9 ~; ?- g* A
    

!using["math","sys"];<font class=\\"jammer\\">0 a6 k/ k\\" V0 z. k! k</font><br /> (:a,b,k,t0)=<br /> <span style=\\"display:none\\">2 s! P; V8 [* w8 s</span> oo{<font class=\\"jammer\\">3 V9 O\\" q7 I. Y2 N6 Q- ?, J& j</font><br />   a=rand[1000,1000], b=rand[1000,1000],<br /> <span style=\\"display:none\\">& t4 i: ^4 g: W6 R</span>   t0=clock(),<br /> <span style=\\"display:none\\">& X4 G% p. k+ t! T</span>   k=a*b,  //矩阵乘<font class=\\"jammer\\">% n+ O9 I$ }4 y+ Z</font><br />   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">0 X8 C# }+ d' J</span> },<font class=\\"jammer\\">7 M# {4 x; D6 C</font><br /> [clock()-t0]/1000;<font class=\\"jammer\\">& l% Q9 ~; ?- g* A</font><br />

结果：

1.         238.447        247.837        247.065        248.105        247.058) z9 |5 P2 {* ]$ T$ z
   
2.         244.123        249.925        247.553        243.981        250.016) d: d. |) g$ U' R( A
   
3.         236.387        252.025        245.651        248.866        248.866
   ; Y8 f8 j+ C& c1 M
4. 2.219 秒
   7 H2 E1 L# l! D) Q0 o& K

复制代码

: Z" [0 s0 L: s" r( k7 Y2、XSLSF（普通的C/C++算法）中的矩阵乘

2. !using["math","sys","XSLSF"];7 B5 m% k, p' b5 t* Z
   
4. (:a,b,k,t0)=
   
5. 2 ^6 |; S3 |# O5 x  {
6. oo{* L' A. M) r: L5 {6 C+ w8 H\\" R
   
8.   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],
   
9. ; {. H$ p  r  l6 s\\" K
10.   t0=clock(),
    
11. 4 `: b  G( R6 ]& l2 @
12.   rmul[k:a,b],  //矩阵乘5 Q+ S1 B8 ~) E. s1 j
    
14.   k[1,3:5,9].outm()
    
15. ) h. X  D; H5 J8 r  N+ w3 }9 @
16. },
    
17. ; f9 p+ J7 V' P, S1 W4 n: |
18. [clock()-t0]/1000;7 A6 J7 D% D\\" C( `: ^4 m) Q& E
    

!using["math","sys","XSLSF"];<font class=\\"jammer\\">7 B5 m% k, p' b5 t* Z</font><br /> (:a,b,k,t0)=<br /> <span style=\\"display:none\\">2 ^6 |; S3 |# O5 x  {</span> oo{<font class=\\"jammer\\">* L' A. M) r: L5 {6 C+ w8 H\\" R</font><br />   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],<br /> <span style=\\"display:none\\">; {. H$ p  r  l6 s\\" K</span>   t0=clock(),<br /> <span style=\\"display:none\\">4 `: b  G( R6 ]& l2 @</span>   rmul[k:a,b],  //矩阵乘<font class=\\"jammer\\">5 Q+ S1 B8 ~) E. s1 j</font><br />   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">) h. X  D; H5 J8 r  N+ w3 }9 @</span> },<br /> <span style=\\"display:none\\">; f9 p+ J7 V' P, S1 W4 n: |</span> [clock()-t0]/1000;<font class=\\"jammer\\">7 A6 J7 D% D\\" C( `: ^4 m) Q& E</font><br />

结果：

1.         262.121        247.583        260.529        259.548        258.328( l0 j, A; m; D/ W/ q! w
   
2.         255.413        246.563        254.356        250.548        251.509
   / P0 u- s5 A4 V. e1 l! F' l2 z
3.         256.152        247.725        259.444        250.827        249.816
   3 d\" \. ~; G0 p) S3 L5 _
4. 10.563 秒; m7 T# x7 P* \& K
   

复制代码