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

( y6 n3 K1 |" y/ M! J- c4 ^

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

9 |; W4 r; D# }9 N1 l% j+ ~

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

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

forcal

例子1代码：
+ Z1 ^2 g$ S0 }3 D; Q( x

2. !using["math"];$ Q# w0 b: t0 A0 o
   
4. mvar:
   
5. & y* J. L, K! K5 q
6. oo{                      //一般在oo函数中调用FcMath函数8 L4 w8 F0 H8 T* R# L& D* l
   
8.   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化# n; e/ X! V3 ]/ d( y. P
   
10.   a.outm(),              //输出矩阵a% k- @2 m( S$ k
    
12.   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出: O3 v( H7 z) z3 _7 l' R4 r
    
14.   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出
    
15. 0 r6 X2 z! H. D' O1 l
16.   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出\\" W. ?; o4 O2 K
    
18. };: g5 O9 D. U' H/ K
    

!using["math"];<font class=\\"jammer\\">$ Q# w0 b: t0 A0 o</font><br /> mvar:<br /> <span style=\\"display:none\\">& y* J. L, K! K5 q</span> oo{                      //一般在oo函数中调用FcMath函数<font class=\\"jammer\\">8 L4 w8 F0 H8 T* R# L& D* l</font><br />   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化<font class=\\"jammer\\"># n; e/ X! V3 ]/ d( y. P</font><br />   a.outm(),              //输出矩阵a<font class=\\"jammer\\">% k- @2 m( S$ k</font><br />   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">: O3 v( H7 z) z3 _7 l' R4 r</font><br />   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出<br /> <span style=\\"display:none\\">0 r6 X2 z! H. D' O1 l</span>   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">\\" W. ?; o4 O2 K</font><br /> };<font class=\\"jammer\\">: g5 O9 D. U' H/ K</font><br />

结果：

1.        0.211319   4.91638e-002       0.144638       0.153259       0.852615
   . C+ h, I2 i1 b
2.        0.630646       0.927048       0.440308       0.162857       0.556854
   ) u2 C' P  }0 O9 q
3.         0.43309        0.34552       0.563919       0.937164       0.209641
   0 r# }( G0 l5 o# N8 r
4.        0.603271       0.727676       0.130951   5.35736e-002       0.197937% d* ?8 A+ x1 `/ z8 z* u
   
5.        0.576004       0.747589   1.17645e-002       0.363892       0.280777
   9 \2 q+ `( F; d  C; x! e* _
6.        0.646454       0.381088        0.58551        0.26387        0.93692* A) e/ b$ K2 S
   
7. ; \/ L( c8 ^; n. F, E! R( M
   
8.        0.1532595 N1 e! f\" E$ F! E  l2 f\" X) h8 M
   
9.        0.162857
   ' ^- }* W# V8 k& S
10.        0.937164
    ( J& Y% U0 k# W9 v  K
11.    5.35736e-002
    : W/ I7 x\" k% c
12.        0.363892\" I6 G7 s( }! B
    
13.         0.26387; @/ g- H7 W( Y7 {# q
    
14. 
    - u. z4 z: W6 }/ s& j
15.        0.603271       0.727676       0.130951   5.35736e-002       0.197937  `0 A/ y' A5 E; B9 T4 ~# e
    
16. 2 k& q& x: j6 q6 e5 u: ~* u* {
    
17.        0.130951   5.35736e-0024 @0 l5 T/ t\" D! ~9 S0 |. ~) X/ `
    
18.    1.17645e-002       0.363892
    2 w  ~1 V* b1 O6 E& R3 a
19.         0.58551        0.26387) |5 x7 B+ L& w' j
    
20. 
    6 S! U# ^1 Y8 D- W( N6 a6 \

复制代码

* ~, d3 K) F% O" ^% p7 Z: G例子2代码：
$ _/ s- c4 s. \- Z2 H
7 J9 w6 z( G. q4 a8 Q* ]  m

2. f(x1,x2,x3,y1,y2,y3)=      //函数定义, U8 q/ C, }- b$ \- j( Z0 U
   
4. {
   
5. / N- l1 r: i& b8 l, ~% H
6.     y1=x1*x1+x2*x2+x3*x3-1.0,
   
7. 7 M* T9 J  ?: Z
8.     y2=2.0*x1*x1+x2*x2-4.0*x3,
   
9.   D9 s9 J5 g- c
10.     y3=3.0*x1*x1-4.0*x2+x3*x36 n* a. `$ t3 ^
    
12. };! K: v8 Y( R3 L1 N) }
    
14. !using["math","sys"];/ z. u: U4 \, C7 V# W0 }
    
16. mvar:
    
17. 4 N* y9 c  x% b: X9 I  x
18. oo{
    
19. 5 p. ]7 r) [, X' K/ R\\" U5 [\\" R8 {
20.   x=array(3),! P8 u% j* ]- l$ k/ C4 b' }
    
22.   x.SA[0 : 1,1,1],       //设置初值为1,1,1# |5 Z, l6 Q; V\\" C
    
24.   i=netn[HFor("f"),x],   //拟牛顿法解方程
    
25. % H4 }- P1 z3 ]
26.   x.outm(),              //输出结果
    
27. 3 i3 o* I; r/ S' C# j' w' B* j
28.   i                      //返回迭代次数
    
29. # T4 k, j3 o& c
30. };
    
31. : J  v, U- v, D' z# P! \4 X

f(x1,x2,x3,y1,y2,y3)=      //函数定义<font class=\\"jammer\\">, U8 q/ C, }- b$ \- j( Z0 U</font><br /> {<br /> <span style=\\"display:none\\">/ N- l1 r: i& b8 l, ~% H</span>     y1=x1*x1+x2*x2+x3*x3-1.0,<br /> <span style=\\"display:none\\">7 M* T9 J  ?: Z</span>     y2=2.0*x1*x1+x2*x2-4.0*x3,<br /> <span style=\\"display:none\\">  D9 s9 J5 g- c</span>     y3=3.0*x1*x1-4.0*x2+x3*x3<font class=\\"jammer\\">6 n* a. `$ t3 ^</font><br /> };<font class=\\"jammer\\">! K: v8 Y( R3 L1 N) }</font><br /> !using["math","sys"];<font class=\\"jammer\\">/ z. u: U4 \, C7 V# W0 }</font><br /> mvar:<br /> <span style=\\"display:none\\">4 N* y9 c  x% b: X9 I  x</span> oo{<br /> <span style=\\"display:none\\">5 p. ]7 r) [, X' K/ R\\" U5 [\\" R8 {</span>   x=array(3),<font class=\\"jammer\\">! P8 u% j* ]- l$ k/ C4 b' }</font><br />   x.SA[0 : 1,1,1],       //设置初值为1,1,1<font class=\\"jammer\\"># |5 Z, l6 Q; V\\" C</font><br />   i=netn[HFor("f"),x],   //拟牛顿法解方程<br /> <span style=\\"display:none\\">% H4 }- P1 z3 ]</span>   x.outm(),              //输出结果<br /> <span style=\\"display:none\\">3 i3 o* I; r/ S' C# j' w' B* j</span>   i                      //返回迭代次数<br /> <span style=\\"display:none\\"># T4 k, j3 o& c</span> };<br /> <span style=\\"display:none\\">: J  v, U- v, D' z# P! \4 X</span>

结果：
  0.785197       0.496611       0.369923

forcal

效率测试：
( Y. N* T; ?3 c6 [& y" A: K+ Jsimwe的网友lin2009 的matlab代码：

1. clear all
     U9 j* o7 t4 c: [4 K3 {
2. clc8 ^3 k6 p2 V. \
   
3. tic
   % }- Q8 g, Q) @9 F1 _7 d) u
4. k = zeros(5,5); % //生成5×5全0矩阵9 C' h# r9 u8 F/ m% w3 }
   
5. % 循环计算以下程序段100000次：
   ( Y\" j2 }( w( v  C
6. for m = 1:100000% C3 I. [\" c6 _7 a
   
7.     a = rand(5,7);
   , g9 ~- {. }3 e5 x( B5 I' U- N
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化/ x6 S# y/ g5 G; U; u( A; ^
   
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   + m) K* v; L' O9 T, ?
10. end
    ; W# I. w) y* L; \9 n
11. k
    . }3 q: b& H/ l7 W5 O% s, X1 ?
12. toc1 ~; U2 [! t8 R, z0 }, S- S$ z
    

复制代码

Forcal代码：
$ w# I4 o  `( y) W! K
运行稍快的代码，比matlab约快10％吧？
& y- `" p: b1 S# A$ O

2. !using["math","sys"];/ _( F* |\\" m2 w/ B  ?\\" U( {
   
4. mvar:2 z6 Y- U, O' o0 ~, @2 a
   
6. t0=clock(),
   
7. , X5 w/ o6 }4 g& ]0 R
8. oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为04 F! E- y7 K9 L8 g- r
   
10. i=0,(i<1000 00).while{ //循环计算1000 00次
    
11. % q+ n  L( U5 a! Z\\" V6 v& G
12.   oo{4 j( E7 i: T' z3 q) _1 `
    
14.     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
    
15. * w1 \* v0 \# H* b1 E# U6 C6 t! p. t
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)) }2 K) B/ m/ Y0 [8 z
    
18.   },, T  Z\\" H6 y- E+ s
    
20.   i++
    
21. . p5 @8 N: U7 U6 g2 W3 M\\" {
22. },4 L; a( o# a/ C% M# ?; k
    
24. k.outm(),             //输出矩阵k，然后销毁k- D' W5 a/ f1 f8 j\\" G
    
26. [clock()-t0]/1000;    //得到计算时间，秒

!using["math","sys"];<font class=\\"jammer\\">/ _( F* |\\" m2 w/ B  ?\\" U( {</font><br /> mvar:<font class=\\"jammer\\">2 z6 Y- U, O' o0 ~, @2 a</font><br /> t0=clock(),<br /> <span style=\\"display:none\\">, X5 w/ o6 }4 g& ]0 R</span> oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0<font class=\\"jammer\\">4 F! E- y7 K9 L8 g- r</font><br /> i=0,(i<1000 00).while{ //循环计算1000 00次<br /> <span style=\\"display:none\\">% q+ n  L( U5 a! Z\\" V6 v& G</span>   oo{<font class=\\"jammer\\">4 j( E7 i: T' z3 q) _1 `</font><br />     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化<br /> <span style=\\"display:none\\">* w1 \* v0 \# H* b1 E# U6 C6 t! p. t</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)<font class=\\"jammer\\">) }2 K) B/ m/ Y0 [8 z</font><br />   },<font class=\\"jammer\\">, T  Z\\" H6 y- E+ s</font><br />   i++<br /> <span style=\\"display:none\\">. p5 @8 N: U7 U6 g2 W3 M\\" {</span> },<font class=\\"jammer\\">4 L; a( o# a/ C% M# ?; k</font><br /> k.outm(),             //输出矩阵k，然后销毁k<font class=\\"jammer\\">- D' W5 a/ f1 f8 j\\" G</font><br /> [clock()-t0]/1000;    //得到计算时间，秒

在我的电脑上运行时间为3.344秒。
$ _6 J: u! M$ z! Q: G" n7 {6 C6 i
- Z4 K& e; `7 i" O9 ?比较好看些的代码，似乎也比matlab稍快吧？
( S; e2 Q2 v& q& Y

2. !using["math","sys"];
   
3. 5 T5 y2 }% c2 Q, N7 }. C, _
4. (:t0,k,i,a,b)=
   
5. 6 G4 Q+ @\\" m$ n2 i* g* |% H/ C\\" x9 R
6. {; F+ e, |; w6 k: k
   
8.   t0=clock(),
   
9. / ?* @3 [: `! y) Y2 |1 x) `/ Q* k
10.   k=zeros[5,5],0 k5 z( A* a, C\\" t0 S
    
12.   i=0,(i<1000 00).while{
    
13. 1 ^% |$ I& t1 h! f( e) ~& {
14.     oo{) w- |! e\\" o% a
    
16.       a=rand[5,7], b=rand[7,5],' B+ ^. c! w6 k
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)' t1 Z# O0 ?: ~- d! t8 T; J
    
20.     },. Z2 U6 S) w3 a, ?- Z
    
22.     i++: m( |\\" Z. O' p& t& {, [
    
24.   },
    
25. - Z\\" G( }+ R/ O1 _6 b\\" J) W
26.   k.outm().delete(),% z; i( T! E1 Z# j6 Z3 b% J
    
28.   [clock()-t0]/1000
    
29. 6 Z\\" n7 T$ o+ h4 D. d6 u% R
30. };

!using["math","sys"];<br /> <span style=\\"display:none\\">5 T5 y2 }% c2 Q, N7 }. C, _</span> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">6 G4 Q+ @\\" m$ n2 i* g* |% H/ C\\" x9 R</span> {<font class=\\"jammer\\">; F+ e, |; w6 k: k</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">/ ?* @3 [: `! y) Y2 |1 x) `/ Q* k</span>   k=zeros[5,5],<font class=\\"jammer\\">0 k5 z( A* a, C\\" t0 S</font><br />   i=0,(i<1000 00).while{<br /> <span style=\\"display:none\\">1 ^% |$ I& t1 h! f( e) ~& {</span>     oo{<font class=\\"jammer\\">) w- |! e\\" o% a</font><br />       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">' B+ ^. c! w6 k</font><br />       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)<font class=\\"jammer\\">' t1 Z# O0 ?: ~- d! t8 T; J</font><br />     },<font class=\\"jammer\\">. Z2 U6 S) w3 a, ?- Z</font><br />     i++<font class=\\"jammer\\">: m( |\\" Z. O' p& t& {, [</font><br />   },<br /> <span style=\\"display:none\\">- Z\\" G( }+ R/ O1 _6 b\\" J) W</span>   k.outm().delete(),<font class=\\"jammer\\">% z; i( T! E1 Z# j6 Z3 b% J</font><br />   [clock()-t0]/1000<br /> <span style=\\"display:none\\">6 Z\\" n7 T$ o+ h4 D. d6 u% R</span> };

2 Q/ f: G5 d5 F3 l在我的电脑上运行时间为3.579秒。
" K* j$ i( c0 ?. q4 T
该例子的理论结果是每个元素均为275000。
+ J3 Z4 r; B+ S2 ]4 W
. g/ j1 a; O! d, N1 ^7 e我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。

forcal

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

2. !using["math"];/ H* w0 }6 N+ t0 s* U3 v3 X
   
4. mvar:! K# V9 s4 w6 \4 {7 G, }
   
6. oo{
   
7. ! l$ `/ r$ @& _; B3 s
8.   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],5 N) q; C+ F\\" k& ?6 S+ p1 E' |1 B
   
10.   a1.outm[5,1,1],
    
11. 4 K\\" E3 D8 u* y' |- E2 O& M7 l
12.   a2.outm[5,1,1],
    
13. & g: d7 Y0 V- S; N) I3 r
14.   a3.outm[5,1,1],3 P9 ]% X5 ], h/ W& t
    
16.   a4.outm[5,1,1],
    
17. & o. ~, B0 C5 P4 |+ e
18.   a=a1+a2+a3+a4,
    
19. 4 L+ N% @/ w. A
20.   a.outm[5,1,1],
    
21. . U3 U: [/ W( j\\" K. @- p4 h
22.   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]0 i' @/ [$ x1 j# a
    
24. };
    
25. # g3 t7 K4 u- e

!using["math"];<font class=\\"jammer\\">/ H* w0 }6 N+ t0 s* U3 v3 X</font><br /> mvar:<font class=\\"jammer\\">! K# V9 s4 w6 \4 {7 G, }</font><br /> oo{<br /> <span style=\\"display:none\\">! l$ `/ r$ @& _; B3 s</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\\">5 N) q; C+ F\\" k& ?6 S+ p1 E' |1 B</font><br />   a1.outm[5,1,1],<br /> <span style=\\"display:none\\">4 K\\" E3 D8 u* y' |- E2 O& M7 l</span>   a2.outm[5,1,1],<br /> <span style=\\"display:none\\">& g: d7 Y0 V- S; N) I3 r</span>   a3.outm[5,1,1],<font class=\\"jammer\\">3 P9 ]% X5 ], h/ W& t</font><br />   a4.outm[5,1,1],<br /> <span style=\\"display:none\\">& o. ~, B0 C5 P4 |+ e</span>   a=a1+a2+a3+a4,<br /> <span style=\\"display:none\\">4 L+ N% @/ w. A</span>   a.outm[5,1,1],<br /> <span style=\\"display:none\\">. U3 U: [/ W( j\\" K. @- p4 h</span>   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]<font class=\\"jammer\\">0 i' @/ [$ x1 j# a</font><br /> };<br /> <span style=\\"display:none\\"># g3 t7 K4 u- e</span>

说明：
linspace(8,12,5)：生成一维数组，共5个元素8～12
& P, `& B3 A- x3 w0 pa1.outm[5,1,1]：输出**数组a1，连下标一起输出
: P6 Z) @5 \# E' [Sum[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)。

2 g9 N, A% }8 a  E+ W; }- P: L结果（最终求和结果是1320）：
2 |8 C6 s) y+ z: e5 M
, f5 X2 j! n+ r(0,0,0,*)              1.            1.            1.            1.            1.
, L; c, u0 w% V1 V% \+ N(0,0,1,*)              1.            1.            1.            1.            1.
7 i+ t$ o2 ^' c( s2 W(0,1,0,*)              1.            1.            1.            1.            1.
) M4 N1 ~! u7 e, y( z8 Z(0,1,1,*)              1.            1.            1.            1.            1.
(0,2,0,*)              1.            1.            1.            1.            1.
(0,2,1,*)              1.            1.            1.            1.            1.
0 r2 f7 Q. y4 ^6 C  l! _& _2 R(1,0,0,*)              2.            2.            2.            2.            2.
(1,0,1,*)              2.            2.            2.            2.            2.
. D$ U9 M% j& U) {7 o& J* L" z(1,1,0,*)              2.            2.            2.            2.            2.
# }! g' Z7 D5 N  w3 ^1 E3 `7 n(1,1,1,*)              2.            2.            2.            2.            2.
! a1 [9 ^! {( E$ _" M+ u(1,2,0,*)              2.            2.            2.            2.            2.
( T: J6 Q# A& G- N9 _( H7 O$ @(1,2,1,*)              2.            2.            2.            2.            2.
+ s/ p9 M% V' b" _3 a) B  W: n* J% y
(0,0,0,*)              3.            3.            3.            3.            3.
" F' x, ]3 ~0 i. L6 l! w( n(0,0,1,*)              3.            3.            3.            3.            3.
" ~! d( \4 ^  V5 U(0,1,0,*)              4.            4.            4.            4.            4.
7 P4 t7 F; i0 f% N# [$ s2 S3 U3 W* e(0,1,1,*)              4.            4.            4.            4.            4.
(0,2,0,*)              5.            5.            5.            5.            5.
0 n& P6 e6 p; h(0,2,1,*)              5.            5.            5.            5.            5.
) N% V# o- q2 X5 h6 r" R(1,0,0,*)              3.            3.            3.            3.            3.
(1,0,1,*)              3.            3.            3.            3.            3.
(1,1,0,*)              4.            4.            4.            4.            4.
(1,1,1,*)              4.            4.            4.            4.            4.
(1,2,0,*)              5.            5.            5.            5.            5.
2 u' S: z- `! a' w) e(1,2,1,*)              5.            5.            5.            5.            5.

6 K2 p1 v, _2 a(0,0,0,*)              6.            6.            6.            6.            6.
# `& M- z# p( B" R  _(0,0,1,*)              7.            7.            7.            7.            7.
7 e, M  S4 u" N# W(0,1,0,*)              6.            6.            6.            6.            6.
7 R3 _0 z, b0 o$ @2 z8 V1 v(0,1,1,*)              7.            7.            7.            7.            7.
) l% g0 {5 m9 y2 I/ h$ e(0,2,0,*)              6.            6.            6.            6.            6.
(0,2,1,*)              7.            7.            7.            7.            7.
(1,0,0,*)              6.            6.            6.            6.            6.
0 u3 d$ d# r) Q(1,0,1,*)              7.            7.            7.            7.            7.
( n$ G, Q! Q( H; C. y3 ~(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.
& C1 I# ?8 H! l3 N(0,0,1,*)              8.            9.           10.           11.           12.
. Y0 E# z" _2 D' C(0,1,0,*)              8.            9.           10.           11.           12.
(0,1,1,*)              8.            9.           10.           11.           12.
(0,2,0,*)              8.            9.           10.           11.           12.
. [- T; r) ?6 D: ^7 A! O. G4 |(0,2,1,*)              8.            9.           10.           11.           12.
9 [$ z8 Z$ k8 [  {; r) L. P(1,0,0,*)              8.            9.           10.           11.           12.
(1,0,1,*)              8.            9.           10.           11.           12.
(1,1,0,*)              8.            9.           10.           11.           12.
(1,1,1,*)              8.            9.           10.           11.           12.
6 Y: I3 o" l$ y8 R(1,2,0,*)              8.            9.           10.           11.           12.
3 |4 W* D# y1 Y# c(1,2,1,*)              8.            9.           10.           11.           12.

$ M( y. S! h2 D* _0 s  ]2 d& z0 z(0,0,0,*)             18.           19.           20.           21.           22.
(0,0,1,*)             19.           20.           21.           22.           23.
) F0 }+ H) H: B( P% s* C(0,1,0,*)             19.           20.           21.           22.           23.
(0,1,1,*)             20.           21.           22.           23.           24.
(0,2,0,*)             20.           21.           22.           23.           24.
(0,2,1,*)             21.           22.           23.           24.           25.
" m4 i, Z# A; f6 S4 L9 \2 X(1,0,0,*)             19.           20.           21.           22.           23.
(1,0,1,*)             20.           21.           22.           23.           24.
$ Y6 v& ~1 k, V; b, V$ y(1,1,0,*)             20.           21.           22.           23.           24.
3 K) o/ H  m" Z, t% j# Z(1,1,1,*)             21.           22.           23.           24.           25.
(1,2,0,*)             21.           22.           23.           24.           25.
8 }9 p0 m. y8 E0 G. u, ]: S1 V5 I7 N2 l(1,2,1,*)             22.           23.           24.           25.           26.
+ N- t) C. p/ }5 p6 J2 o
' t, A6 T3 z4 |1 P/ E2 m(0,0,*)            100.          105.
- J+ p+ o+ C7 `- g" f, x(0,1,*)            105.          110.
(0,2,*)            110.          115.
(1,0,*)            105.          110.
(1,1,*)            110.          115.
(1,2,*)            115.          120.

; F5 j. d7 R4 z- {$ o/ g(0,*)            205.          215.          225.
% m& N1 _4 h& p9 M* m" u(1,*)            215.          225.          235.
. o/ Z) D, P+ `( i9 c# }# T
(0,*)            645.
(1,*)            675.
  n3 Y; i$ {' ~) [1 e
1320.

forcal

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

2. !using["math","sys"];$ D1 B\\" z; c/ Y$ l4 f
   
4. mvar:
   
5. $ x# E; B7 x6 s+ d% D5 r+ ^
6. (:p1,p2,p3,a,b)=  |4 Z, ~\\" O  @' R
   
8. {
   
9. ) t( x$ U0 t* n4 w7 r$ \* q! `( x
10.   oo{
    
11. 1 p* A- M; D' t
12.     a=array[1000].rand(),! F, }6 Q7 ?+ x0 q! B; \
    
14.     b=array[1000].rand(),
    
15. * q* b1 z- |6 E9 m8 s  f
16.     p1=array[1000,1000],) x: `$ N. F\\" j* D. C# U
    
18.     p2=array[1000,1000],8 d  q) G: u% J: _0 a\\" y. d2 A
    
20.     p3=array[1000,1000],
    
21. ) n9 g% R. L) _2 n
22.     t0=clock(),9 y7 {  ~. e+ b$ h
    
24.     ndgrid(a,b,&A,&B),
    
25. 6 p2 c# G9 ?( d: X- b' l* @
26.     p1.=A+B$ f6 x& A8 ~7 X
    
28.   },
    
29. 4 q5 Q9 e4 X$ X. P3 t# W/ m# S. g; S
30.   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},: y( z+ |& ^; [; f, l
    
32.   lena=FCDLen(a),
    
33. 6 s0 @' h  \7 t: h' ^# F\\" Q5 c
34.   lenb=FCDLen(b),
    
35. - {+ ]7 K: q( m' C
36.   t0=clock(),$ ~- U  h1 @% e/ b0 q3 M% j\\" k
    
38.   m = lenb-1, (m>=0).while{9 B4 |8 X: Q$ B0 t: N
    
40.     oo{p2(m,neg) = a+rn[b(m)]},
    
41. % I( O  h* w7 o' _! B
42.     m--# ]( O5 m5 D2 n\\" O
    
44.   },
    
45. 1 }! o! K8 D  p9 l% N9 H: A, ^3 t- {
46.   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},
    
47. ; K9 x# R* A: m4 e. D% C! z1 \
48.   t0=clock(),/ Q\\" r1 N9 F/ R+ o$ K( S3 O
    
50.   m = lenb-1, (m>=0).while{: \# h: U2 t% H. q& P; f
    
52.     n = lena-1, (n>=0).while{$ ]- ^0 M9 {\\" N1 z3 `, p
    
54.       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些8 N: ~) O6 B2 N. r# a
    
56.       A(p3,m,n) = A(a,n)+A(b,m),\\" y0 Y) f$ ^- Y; {9 {
    
58.       n--4 l- i! s0 l$ `- N9 h
    
60.     },# D; W3 w( N' D\\" p- h( S+ L
    
62.     m--
    
63. + T. b0 L9 L) n. {
64.   },6 G/ B4 t% d( E6 A* p' R+ F: |
    
66.   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
    
67. & O. L6 z$ I6 V
68. };; i9 H. N5 r8 c- L# g& t( G
    

!using["math","sys"];<font class=\\"jammer\\">$ D1 B\\" z; c/ Y$ l4 f</font><br /> mvar:<br /> <span style=\\"display:none\\">$ x# E; B7 x6 s+ d% D5 r+ ^</span> (:p1,p2,p3,a,b)=<font class=\\"jammer\\">  |4 Z, ~\\" O  @' R</font><br /> {<br /> <span style=\\"display:none\\">) t( x$ U0 t* n4 w7 r$ \* q! `( x</span>   oo{<br /> <span style=\\"display:none\\">1 p* A- M; D' t</span>     a=array[1000].rand(),<font class=\\"jammer\\">! F, }6 Q7 ?+ x0 q! B; \</font><br />     b=array[1000].rand(),<br /> <span style=\\"display:none\\">* q* b1 z- |6 E9 m8 s  f</span>     p1=array[1000,1000],<font class=\\"jammer\\">) x: `$ N. F\\" j* D. C# U</font><br />     p2=array[1000,1000],<font class=\\"jammer\\">8 d  q) G: u% J: _0 a\\" y. d2 A</font><br />     p3=array[1000,1000],<br /> <span style=\\"display:none\\">) n9 g% R. L) _2 n</span>     t0=clock(),<font class=\\"jammer\\">9 y7 {  ~. e+ b$ h</font><br />     ndgrid(a,b,&A,&B),<br /> <span style=\\"display:none\\">6 p2 c# G9 ?( d: X- b' l* @</span>     p1.=A+B<font class=\\"jammer\\">$ f6 x& A8 ~7 X</font><br />   },<br /> <span style=\\"display:none\\">4 q5 Q9 e4 X$ X. P3 t# W/ m# S. g; S</span>   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},<font class=\\"jammer\\">: y( z+ |& ^; [; f, l</font><br />   lena=FCDLen(a),<br /> <span style=\\"display:none\\">6 s0 @' h  \7 t: h' ^# F\\" Q5 c</span>   lenb=FCDLen(b),<br /> <span style=\\"display:none\\">- {+ ]7 K: q( m' C</span>   t0=clock(),<font class=\\"jammer\\">$ ~- U  h1 @% e/ b0 q3 M% j\\" k</font><br />   m = lenb-1, (m>=0).while{<font class=\\"jammer\\">9 B4 |8 X: Q$ B0 t: N</font><br />     oo{p2(m,neg) = a+rn[b(m)]},<br /> <span style=\\"display:none\\">% I( O  h* w7 o' _! B</span>     m--<font class=\\"jammer\\"># ]( O5 m5 D2 n\\" O</font><br />   },<br /> <span style=\\"display:none\\">1 }! o! K8 D  p9 l% N9 H: A, ^3 t- {</span>   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},<br /> <span style=\\"display:none\\">; K9 x# R* A: m4 e. D% C! z1 \</span>   t0=clock(),<font class=\\"jammer\\">/ Q\\" r1 N9 F/ R+ o$ K( S3 O</font><br />   m = lenb-1, (m>=0).while{<font class=\\"jammer\\">: \# h: U2 t% H. q& P; f</font><br />     n = lena-1, (n>=0).while{<font class=\\"jammer\\">$ ]- ^0 M9 {\\" N1 z3 `, p</font><br />       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些<font class=\\"jammer\\">8 N: ~) O6 B2 N. r# a</font><br />       A(p3,m,n) = A(a,n)+A(b,m),<font class=\\"jammer\\">\\" y0 Y) f$ ^- Y; {9 {</font><br />       n--<font class=\\"jammer\\">4 l- i! s0 l$ `- N9 h</font><br />     },<font class=\\"jammer\\"># D; W3 w( N' D\\" p- h( S+ L</font><br />     m--<br /> <span style=\\"display:none\\">+ T. b0 L9 L) n. {</span>   },<font class=\\"jammer\\">6 G/ B4 t% d( E6 A* p' R+ F: |</font><br />   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}<br /> <span style=\\"display:none\\">& O. L6 z$ I6 V</span> };<font class=\\"jammer\\">; i9 H. N5 r8 c- L# g& t( G</font><br />

结果：
ndgrid: 3.2001e-002
$ t5 z; J9 D6 U  U  ], r2 v7 ]6 Efor1: 1.4999e-002
for2: 1.86
! R9 \( T; A8 d$ f' H) a0 ]" @. `) F5 T
6 w6 n) n4 Z$ h& d( X* A

forcal

一段程序的Forcal实现：
" `; B$ o* ^* X% G% U3 J
9 b% N# Q* W: S! T& l& i//用C++代码描述为：
3 |' h# ^% [0 O) m) w8 Vs=0.0;
for(x=0.0;x<=1.0;x=x+0.0011)
5 ]% W- I8 A# S# E6 L) @{
, N5 _7 M  [8 P   for(y=1.0;y<=2.0;y=y+0.0009)
   {
" |% U. s# ]- d     s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
- i" }) f" f, J* [   }
/ P9 r. H" |9 h  T1 X5 x2 T# N}
, V) Y" y8 T7 l& L3 l
" O3 O3 d- _  {1、**数组求和函数Sum

2. !using["math","sys"];: D; {4 T1 y  A: t1 A8 Z$ S
   
4. mvar:) V  |9 _8 D+ J% ]
   
6. t=clock(),
   
7. 2 C/ b+ w4 i+ e: c9 u
8. oo{
   
9. 3 N' o. f* |1 q
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],/ n3 t$ J1 u; b6 u( Z$ h. K3 \& u2 ~
    
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. \\" y# }' }) S0 ~, L: j& _% |
14. };3 t& d) u2 D. [
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">: D; {4 T1 y  A: t1 A8 Z$ S</font><br /> mvar:<font class=\\"jammer\\">) V  |9 _8 D+ J% ]</font><br /> t=clock(),<br /> <span style=\\"display:none\\">2 C/ b+ w4 i+ e: c9 u</span> oo{<br /> <span style=\\"display:none\\">3 N' o. f* |1 q</span>   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<font class=\\"jammer\\">/ n3 t$ J1 u; b6 u( Z$ h. K3 \& u2 ~</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\\">\\" y# }' }) S0 ~, L: j& _% |</span> };<font class=\\"jammer\\">3 t& d) u2 D. [</font><br /> [clock()-t]/1000;

结果：
1008606.64947441
9 I# ~5 A0 Z: I2 }0.625   //时间
. i/ F9 |  p9 `9 Z8 M* y
2、求和函数sum

  R) `3 [" a( u! n9 Z# |7 {  \

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   ! g\" i) _: Z5 p6 H! M% t! [
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

) j  l  j9 R0 v- l6 o! W* j结果：
1008606.64947441
0.719   //时间

3 T6 ]: q+ q# a9 P% ]! L3、while循环

1. mvar:, `, N9 g' W* |/ F! x
   
2. t=sys::clock();
   1 x: J9 `: Q7 d$ s  o
3. s=0,x=0, \" y2 i  @7 m% M! g! z\" A3 n
   
4. while{x<=1,  //while循环算法； # h) Z: K, z, v$ {' \; J
   
5.    y=1,
   6 B0 q* |& z$ m3 B8 B0 F
6.    while{y<=2,
   \" i/ k4 }3 f8 y1 ~& [3 n% g# F
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), & C6 \! n6 U: G& S$ J: v
   
8.        y=y+0.0009 ' E% W& C  ?# N' \1 ?$ W
   
9.       },
   + D  {! F5 `2 A7 ?5 U\" J
10.    x=x+0.0011 7 N\" Y8 |: E& e
    
11. },
    + q) E4 O. R) S. q5 x2 H+ q* Y
12. s;
    . h5 w! _+ G+ Q- l
13. [sys::clock()-t]/1000;

复制代码

  j/ I; P% P  L% p9 s% s0 b结果：
) f/ ~) c4 L+ j: q2 E+ J$ |1008606.64947441
1 v/ F+ M, @6 ?0.734   //时间
# i1 L  Y- ~# E% |

qbist

好深奥！~~~~

forcal

好深奥！~~~~
: x9 }' Q9 }! [' ]qbist 发表于 2010-10-7 14:56

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

qbist

回复 forcal 的帖子

$ N5 N( f- }0 S
% e9 y) [  t; w4 O& u  g; I5 j! D1 [    嗯！！！

forcal

改进了FcMath中的矩阵乘算法，不知与matlab还有多大差距，朋友们可帮助测一下。
; D8 I9 l% n1 w, {9 Z
7 q5 [* y$ f! [  X以下是FcMath中的矩阵乘与徐士良算法库XSLSF（普通的C/C++算法）中的矩阵乘的效率比较：
. p) V0 i& m' X6 J3 t+ I1 Z1 o
1、FcMath中的矩阵乘
# P3 S# m$ i% t9 Z! L

2. !using["math","sys"];/ U' e\\" I2 a' Y2 {0 x
   
4. (:a,b,k,t0)=
   
5. . p* A! g. r2 o3 W, p0 Z\\" y
6. oo{0 f! F' ^  {( ^6 `: G
   
8.   a=rand[1000,1000], b=rand[1000,1000],
   
9. $ o( \- u8 x$ q# m9 N0 r
10.   t0=clock(),
    
11. 2 d' T5 \! Z% e
12.   k=a*b,  //矩阵乘
    
13. ! h% ]! \5 E& B2 ^
14.   k[1,3:5,9].outm()
    
15. - l+ @& O: c+ t
16. },
    
17. 3 g& }: d8 y. d4 o9 |1 |
18. [clock()-t0]/1000;! ?: @* [3 |/ Y' u! C+ m% G$ d8 g
    

!using["math","sys"];<font class=\\"jammer\\">/ U' e\\" I2 a' Y2 {0 x</font><br /> (:a,b,k,t0)=<br /> <span style=\\"display:none\\">. p* A! g. r2 o3 W, p0 Z\\" y</span> oo{<font class=\\"jammer\\">0 f! F' ^  {( ^6 `: G</font><br />   a=rand[1000,1000], b=rand[1000,1000],<br /> <span style=\\"display:none\\">$ o( \- u8 x$ q# m9 N0 r</span>   t0=clock(),<br /> <span style=\\"display:none\\">2 d' T5 \! Z% e</span>   k=a*b,  //矩阵乘<br /> <span style=\\"display:none\\">! h% ]! \5 E& B2 ^</span>   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">- l+ @& O: c+ t</span> },<br /> <span style=\\"display:none\\">3 g& }: d8 y. d4 o9 |1 |</span> [clock()-t0]/1000;<font class=\\"jammer\\">! ?: @* [3 |/ Y' u! C+ m% G$ d8 g</font><br />

结果：
! U# |" o. Z2 z5 v- F3 T

1.         238.447        247.837        247.065        248.105        247.0586 e& k6 \3 ~* U' F6 `
   
2.         244.123        249.925        247.553        243.981        250.016% T' d+ _* U; |' b8 D/ N
   
3.         236.387        252.025        245.651        248.866        248.8666 p* O' }5 w2 u& A( ~
   
4. 2.219 秒
   0 W9 A$ w: g0 ?% Z0 B' ?0 c) W6 I6 u' P

复制代码

2、XSLSF（普通的C/C++算法）中的矩阵乘

2. !using["math","sys","XSLSF"];1 w( t$ X2 H4 X5 ~0 t
   
4. (:a,b,k,t0)=
   
5. , Y4 d$ q& _% ~\\" ^9 a% l4 B( v
6. oo{% W2 R; Z& E6 f2 P
   
8.   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],/ I8 W% L2 {+ W7 _
   
10.   t0=clock(),
    
11. # T+ b0 c5 E! z3 b, n+ N8 e; Y* {. \
12.   rmul[k:a,b],  //矩阵乘
    
13. / m8 Q; ~$ G' W6 W: \' `
14.   k[1,3:5,9].outm()
    
15. . n  o+ O5 H7 o4 c4 f+ P5 n; [- I2 B
16. },
    
17. + C3 }% q! c& o# {+ v
18. [clock()-t0]/1000;+ n0 U0 Q$ {! G
    

!using["math","sys","XSLSF"];<font class=\\"jammer\\">1 w( t$ X2 H4 X5 ~0 t</font><br /> (:a,b,k,t0)=<br /> <span style=\\"display:none\\">, Y4 d$ q& _% ~\\" ^9 a% l4 B( v</span> oo{<font class=\\"jammer\\">% W2 R; Z& E6 f2 P</font><br />   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],<font class=\\"jammer\\">/ I8 W% L2 {+ W7 _</font><br />   t0=clock(),<br /> <span style=\\"display:none\\"># T+ b0 c5 E! z3 b, n+ N8 e; Y* {. \</span>   rmul[k:a,b],  //矩阵乘<br /> <span style=\\"display:none\\">/ m8 Q; ~$ G' W6 W: \' `</span>   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">. n  o+ O5 H7 o4 c4 f+ P5 n; [- I2 B</span> },<br /> <span style=\\"display:none\\">+ C3 }% q! c& o# {+ v</span> [clock()-t0]/1000;<font class=\\"jammer\\">+ n0 U0 Q$ {! G</font><br />

结果：

1.         262.121        247.583        260.529        259.548        258.328
   ' P* N' Q/ c1 Z% X5 j
2.         255.413        246.563        254.356        250.548        251.509
   . d. U+ A# A! F1 J
3.         256.152        247.725        259.444        250.827        249.816
   8 a) u0 T/ J- b. |: d8 |
4. 10.563 秒8 g& T, Z8 B: F& ?
   

复制代码

$ `- X7 z9 ?; S6 h$ P* m$ Q. M' B9 [$ A