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

# F8 P3 y! K1 N8 z

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

7 e' W; e3 l. v2 t( x- I! @

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

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

forcal

例子1代码：
( \. y6 P1 F* m2 x/ O7 C

2. !using["math"];4 c) Q) ?+ |2 \: i+ F
   
4. mvar:
   
5. + @4 I) R8 V+ H\\" i
6. oo{                      //一般在oo函数中调用FcMath函数
   
7. 9 H, r9 ]1 J' n+ a% g5 X( `
8.   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化
   
9. 7 O% R; W, c1 Z6 `# X% G) L& G
10.   a.outm(),              //输出矩阵a
    
11. ) D  B$ `/ F! ]$ b9 n7 }6 t! s
12.   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出
    
13. ( k' R: }6 ?1 T( _& C6 |& s
14.   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出9 Z  O0 G9 o+ o0 E
    
16.   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出) |6 |9 H4 q) f  @- X
    
18. };: s5 ]/ v0 Q' m& H: G7 R
    

!using["math"];<font class=\\"jammer\\">4 c) Q) ?+ |2 \: i+ F</font><br /> mvar:<br /> <span style=\\"display:none\\">+ @4 I) R8 V+ H\\" i</span> oo{                      //一般在oo函数中调用FcMath函数<br /> <span style=\\"display:none\\">9 H, r9 ]1 J' n+ a% g5 X( `</span>   a=rand[6,5],           //生成6×5矩阵a，用0～1之间随机数初始化<br /> <span style=\\"display:none\\">7 O% R; W, c1 Z6 `# X% G) L& G</span>   a.outm(),              //输出矩阵a<br /> <span style=\\"display:none\\">) D  B$ `/ F! ]$ b9 n7 }6 t! s</span>   a.subg(neg:3).outm(),  //取矩阵a第4列所有元素组成子矩阵，并输出<br /> <span style=\\"display:none\\">( k' R: }6 ?1 T( _& C6 |& s</span>   a.subg(3:neg).outm(),  //取矩阵a第4行所有元素组成子矩阵，并输出<font class=\\"jammer\\">9 Z  O0 G9 o+ o0 E</font><br />   a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出<font class=\\"jammer\\">) |6 |9 H4 q) f  @- X</font><br /> };<font class=\\"jammer\\">: s5 ]/ v0 Q' m& H: G7 R</font><br />

结果：
- f: V1 ?, I6 H7 u5 E9 V

1.        0.211319   4.91638e-002       0.144638       0.153259       0.852615
   6 d7 u\" n2 F/ t( c& X, k\" n! w* H
2.        0.630646       0.927048       0.440308       0.162857       0.556854
   . T- X4 I+ n/ w6 {) r
3.         0.43309        0.34552       0.563919       0.937164       0.209641\" p- q- e7 ?1 [0 g0 o5 d, P
   
4.        0.603271       0.727676       0.130951   5.35736e-002       0.197937( F5 m# \* ^3 m
   
5.        0.576004       0.747589   1.17645e-002       0.363892       0.280777* B6 O7 Z# `/ k1 I4 F% L3 G7 Q
   
6.        0.646454       0.381088        0.58551        0.26387        0.936922 x: Y/ N$ m$ H4 n
   
7. 
   & }# b% }: R' o# A
8.        0.1532593 W6 M: n9 ?8 W7 @8 E) F
   
9.        0.162857\" B5 B6 _/ P1 J$ C: j; U% A
   
10.        0.937164
    ) r. d, ?\" A2 G  g
11.    5.35736e-002
    5 t8 |7 ?1 F( _/ j6 s. [6 z
12.        0.363892' T6 B7 y. b* |1 q
    
13.         0.26387
    \" J/ c. Y) @6 G# d3 Z: c3 Y
14. 8 n$ i2 }- i\" \; M9 z- H9 I# u
    
15.        0.603271       0.727676       0.130951   5.35736e-002       0.197937
    / c7 v4 T* o. D1 U' b, x0 s
16. 
    + K, _6 T! y( Q$ N. r5 v
17.        0.130951   5.35736e-002
    ' ]! G7 V' P+ p2 p
18.    1.17645e-002       0.363892+ M- e6 d* j9 o, X) [9 w- o/ \
    
19.         0.58551        0.26387
    8 o- U9 R  p5 q2 `\" E
20. 
    1 L0 h+ C) C3 [5 C1 A# ]

复制代码

+ B4 J5 ~! C* H例子2代码：

+ M6 `4 E$ Z" X9 i+ _

2. f(x1,x2,x3,y1,y2,y3)=      //函数定义
   
3. : a& f% U: Y( A' X
4. {4 l- a; ^' M4 T- |6 b+ l6 \
   
6.     y1=x1*x1+x2*x2+x3*x3-1.0,
   
7. 7 L$ f- k4 u7 r! F# P6 N
8.     y2=2.0*x1*x1+x2*x2-4.0*x3,+ q7 y' `/ l' R( A( {
   
10.     y3=3.0*x1*x1-4.0*x2+x3*x3! ]/ ?$ v% R, k3 j9 I9 B- n
    
12. };
    
13. 3 a0 K- v. s0 i; Y
14. !using["math","sys"];
    
15. $ v( ]  u2 j# s
16. mvar:
    
17. : i) L: l/ G1 s
18. oo{3 o\\" s8 y, l7 f5 r
    
20.   x=array(3),
    
21. 3 E% j# m' d6 q1 m1 i; V: F$ u
22.   x.SA[0 : 1,1,1],       //设置初值为1,1,11 Q* B' X; B: L, s# V\\" ?1 ]6 @
    
24.   i=netn[HFor("f"),x],   //拟牛顿法解方程$ O7 F2 a4 v- g1 T( n
    
26.   x.outm(),              //输出结果
    
27. 4 X9 t, T9 t: Y& Z! J
28.   i                      //返回迭代次数. j: d0 w2 j+ N! c9 i
    
30. };
    
31. , W0 q9 I2 {: H9 a

f(x1,x2,x3,y1,y2,y3)=      //函数定义<br /> <span style=\\"display:none\\">: a& f% U: Y( A' X</span> {<font class=\\"jammer\\">4 l- a; ^' M4 T- |6 b+ l6 \</font><br />     y1=x1*x1+x2*x2+x3*x3-1.0,<br /> <span style=\\"display:none\\">7 L$ f- k4 u7 r! F# P6 N</span>     y2=2.0*x1*x1+x2*x2-4.0*x3,<font class=\\"jammer\\">+ q7 y' `/ l' R( A( {</font><br />     y3=3.0*x1*x1-4.0*x2+x3*x3<font class=\\"jammer\\">! ]/ ?$ v% R, k3 j9 I9 B- n</font><br /> };<br /> <span style=\\"display:none\\">3 a0 K- v. s0 i; Y</span> !using["math","sys"];<br /> <span style=\\"display:none\\">$ v( ]  u2 j# s</span> mvar:<br /> <span style=\\"display:none\\">: i) L: l/ G1 s</span> oo{<font class=\\"jammer\\">3 o\\" s8 y, l7 f5 r</font><br />   x=array(3),<br /> <span style=\\"display:none\\">3 E% j# m' d6 q1 m1 i; V: F$ u</span>   x.SA[0 : 1,1,1],       //设置初值为1,1,1<font class=\\"jammer\\">1 Q* B' X; B: L, s# V\\" ?1 ]6 @</font><br />   i=netn[HFor("f"),x],   //拟牛顿法解方程<font class=\\"jammer\\">$ O7 F2 a4 v- g1 T( n</font><br />   x.outm(),              //输出结果<br /> <span style=\\"display:none\\">4 X9 t, T9 t: Y& Z! J</span>   i                      //返回迭代次数<font class=\\"jammer\\">. j: d0 w2 j+ N! c9 i</font><br /> };<br /> <span style=\\"display:none\\">, W0 q9 I2 {: H9 a</span>

结果：
  0.785197       0.496611       0.369923
  }  L2 f  z$ }

forcal

效率测试：
  N; T4 `0 D+ u# `simwe的网友lin2009 的matlab代码：

1. clear all, j, v\" B/ V3 x\" A  e% k% X
   
2. clc\" o2 i% y$ \# p3 x7 O# s7 D4 Z
   
3. tic
   \" C( W) |: u2 m\" `- ^
4. k = zeros(5,5); % //生成5×5全0矩阵+ J% ^) U. u. P- O; C+ K1 H2 P
   
5. % 循环计算以下程序段100000次：6 `. H( I# h( E. M
   
6. for m = 1:1000003 |( |6 N9 l\" n
   
7.     a = rand(5,7);
   3 f* E1 m; A. ^/ Z
8.     b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
   6 x- ]7 g* z' |/ w, S' g
9.     k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
   9 }! x9 w/ f\" V& f% ~
10. end. q5 X0 F/ N- D/ X0 v; T
    
11. k
    , U$ M/ c  n. P
12. toc% m5 a2 s% F0 v* q
    

复制代码

2 H) f% n9 |$ K5 PForcal代码：
' I% I2 Q) j( m* T5 W4 p( Y* |. I
运行稍快的代码，比matlab约快10％吧？

2. !using["math","sys"];
   
3. 6 e( ?9 |& R/ e; y% C
4. mvar:. E0 D; x$ O2 z
   
6. t0=clock(),
   
7. / q! P& c% p2 n# s
8. oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0
   
9. 8 T5 m* _1 K; d' x- B
10. i=0,(i<1000 00).while{ //循环计算1000 00次
    
11. \\" H1 W0 S+ _/ n
12.   oo{/ L5 Z+ @$ P( M
    
14.     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
    
15. 1 R) H6 f1 [: C: r. ~: v, K
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)5 K! d2 R9 n* C: N
    
18.   },6 @' _! A! H# W) K
    
20.   i++. t# s3 p' i; C5 r; e
    
22. },
    
23. $ ?# ~% x) _# T
24. k.outm(),             //输出矩阵k，然后销毁k5 f, k# T9 h) e
    
26. [clock()-t0]/1000;    //得到计算时间，秒

!using["math","sys"];<br /> <span style=\\"display:none\\">6 e( ?9 |& R/ e; y% C</span> mvar:<font class=\\"jammer\\">. E0 D; x$ O2 z</font><br /> t0=clock(),<br /> <span style=\\"display:none\\">/ q! P& c% p2 n# s</span> oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0<br /> <span style=\\"display:none\\">8 T5 m* _1 K; d' x- B</span> i=0,(i<1000 00).while{ //循环计算1000 00次<br /> <span style=\\"display:none\\">\\" H1 W0 S+ _/ n</span>   oo{<font class=\\"jammer\\">/ L5 Z+ @$ P( M</font><br />     a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化<br /> <span style=\\"display:none\\">1 R) H6 f1 [: C: r. ~: v, K</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\\">5 K! d2 R9 n* C: N</font><br />   },<font class=\\"jammer\\">6 @' _! A! H# W) K</font><br />   i++<font class=\\"jammer\\">. t# s3 p' i; C5 r; e</font><br /> },<br /> <span style=\\"display:none\\">$ ?# ~% x) _# T</span> k.outm(),             //输出矩阵k，然后销毁k<font class=\\"jammer\\">5 f, k# T9 h) e</font><br /> [clock()-t0]/1000;    //得到计算时间，秒

0 T; h8 N- Y6 a8 B在我的电脑上运行时间为3.344秒。
) W8 G9 X% b3 b
比较好看些的代码，似乎也比matlab稍快吧？

2. !using["math","sys"];5 e6 e% I0 V/ u* Z* ^\\" B& d
   
4. (:t0,k,i,a,b)=
   
5. \\" z. m) g9 B* U& z) J9 C( K  J
6. {/ k7 Y0 O9 n- v0 R5 L7 t) S0 s
   
8.   t0=clock(),
   
9. ( R6 w* t% U$ M( g4 B! x
10.   k=zeros[5,5],
    
11. 3 T3 M4 ?9 p% r0 Z) T8 w
12.   i=0,(i<1000 00).while{$ d' O8 `/ [- n$ H$ G
    
14.     oo{. B! h9 y- N# P' I
    
16.       a=rand[5,7], b=rand[7,5],) v7 w* `- A3 Y
    
18.       k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
    
19. * ?9 R2 e4 x# R6 K3 ]
20.     },& F8 l7 z3 @0 W. v2 P6 W. N
    
22.     i++
    
23. % r' M, c  K\\" L- B% e( q8 r! b/ [* L
24.   },
    
25. 9 O9 f# B4 {) {% ]) e
26.   k.outm().delete(),+ ?$ i/ B/ S/ r6 F8 ?  [
    
28.   [clock()-t0]/1000
    
29. \\" @7 j1 {1 {1 D/ W! Z1 Q
30. };

!using["math","sys"];<font class=\\"jammer\\">5 e6 e% I0 V/ u* Z* ^\\" B& d</font><br /> (:t0,k,i,a,b)=<br /> <span style=\\"display:none\\">\\" z. m) g9 B* U& z) J9 C( K  J</span> {<font class=\\"jammer\\">/ k7 Y0 O9 n- v0 R5 L7 t) S0 s</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">( R6 w* t% U$ M( g4 B! x</span>   k=zeros[5,5],<br /> <span style=\\"display:none\\">3 T3 M4 ?9 p% r0 Z) T8 w</span>   i=0,(i<1000 00).while{<font class=\\"jammer\\">$ d' O8 `/ [- n$ H$ G</font><br />     oo{<font class=\\"jammer\\">. B! h9 y- N# P' I</font><br />       a=rand[5,7], b=rand[7,5],<font class=\\"jammer\\">) v7 w* `- A3 Y</font><br />       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\\">* ?9 R2 e4 x# R6 K3 ]</span>     },<font class=\\"jammer\\">& F8 l7 z3 @0 W. v2 P6 W. N</font><br />     i++<br /> <span style=\\"display:none\\">% r' M, c  K\\" L- B% e( q8 r! b/ [* L</span>   },<br /> <span style=\\"display:none\\">9 O9 f# B4 {) {% ]) e</span>   k.outm().delete(),<font class=\\"jammer\\">+ ?$ i/ B/ S/ r6 F8 ?  [</font><br />   [clock()-t0]/1000<br /> <span style=\\"display:none\\">\\" @7 j1 {1 {1 D/ W! Z1 Q</span> };

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

" p0 }! R" N7 n( T; `# }该例子的理论结果是每个元素均为275000。

3 l0 f  U( _6 N  R/ ^/ w8 F9 K( |9 C我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。
2 Z; v& F/ k/ M- Q- |& c5 r* J

forcal

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

2. !using["math"];' s! x1 N9 \/ i. M
   
4. mvar:! L3 v4 H9 H\\" o/ O8 o/ X' Y
   
6. oo{. b: A+ x* Z2 S; K& n+ z
   
8.   ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
   
9. * c, d7 F3 S9 V* l
10.   a1.outm[5,1,1],6 G% U. g/ f) s. B0 d& T
    
12.   a2.outm[5,1,1],% v! Q2 a+ b: i; L/ [2 v
    
14.   a3.outm[5,1,1],! Q% G\\" s, s' d/ X9 }# f
    
16.   a4.outm[5,1,1],  O1 N8 h, d7 _3 p) \\\" ?! k
    
18.   a=a1+a2+a3+a4,
    
19. 1 B. ~% R( A5 |3 x! I1 C4 A
20.   a.outm[5,1,1],
    
21. * q  p4 A0 P2 D' p% ]
22.   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]
    
23. 7 S( T8 Y. N, g
24. };
    
25. ! X# S& n& E3 \/ h. e* L

!using["math"];<font class=\\"jammer\\">' s! x1 N9 \/ i. M</font><br /> mvar:<font class=\\"jammer\\">! L3 v4 H9 H\\" o/ O8 o/ X' Y</font><br /> oo{<font class=\\"jammer\\">. b: A+ x* Z2 S; K& n+ z</font><br />   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\\">* c, d7 F3 S9 V* l</span>   a1.outm[5,1,1],<font class=\\"jammer\\">6 G% U. g/ f) s. B0 d& T</font><br />   a2.outm[5,1,1],<font class=\\"jammer\\">% v! Q2 a+ b: i; L/ [2 v</font><br />   a3.outm[5,1,1],<font class=\\"jammer\\">! Q% G\\" s, s' d/ X9 }# f</font><br />   a4.outm[5,1,1],<font class=\\"jammer\\">  O1 N8 h, d7 _3 p) \\\" ?! k</font><br />   a=a1+a2+a3+a4,<br /> <span style=\\"display:none\\">1 B. ~% R( A5 |3 x! I1 C4 A</span>   a.outm[5,1,1],<br /> <span style=\\"display:none\\">* q  p4 A0 P2 D' p% ]</span>   Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]<br /> <span style=\\"display:none\\">7 S( T8 Y. N, g</span> };<br /> <span style=\\"display:none\\">! X# S& n& E3 \/ h. e* L</span>

5 E! I; i) t3 J4 p+ i5 o  R+ i说明：
linspace(8,12,5)：生成一维数组，共5个元素8～12
- z4 p- Q* m' `' x, \' m/ Ua1.outm[5,1,1]：输出**数组a1，连下标一起输出
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)。

结果（最终求和结果是1320）：
/ H! h' X" |; f* }3 U5 [8 J
% e: h' S1 b) c4 O9 O(0,0,0,*)              1.            1.            1.            1.            1.
(0,0,1,*)              1.            1.            1.            1.            1.
1 G: v1 |& i, J5 x; u(0,1,0,*)              1.            1.            1.            1.            1.
(0,1,1,*)              1.            1.            1.            1.            1.
(0,2,0,*)              1.            1.            1.            1.            1.
( U. s  u; [7 f2 W0 [7 {(0,2,1,*)              1.            1.            1.            1.            1.
(1,0,0,*)              2.            2.            2.            2.            2.
" L- M4 \: J& k& h8 e3 h7 R; O8 F7 @/ W(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 l/ ^: m4 ?, C) z(1,2,0,*)              2.            2.            2.            2.            2.
(1,2,1,*)              2.            2.            2.            2.            2.

. |  O7 B& c" J6 F! r(0,0,0,*)              3.            3.            3.            3.            3.
9 V! ]! T% i7 L, ]0 V(0,0,1,*)              3.            3.            3.            3.            3.
(0,1,0,*)              4.            4.            4.            4.            4.
(0,1,1,*)              4.            4.            4.            4.            4.
2 t2 i  g/ l* G( _2 q(0,2,0,*)              5.            5.            5.            5.            5.
: P3 n6 C3 p% }& K0 X( z6 |(0,2,1,*)              5.            5.            5.            5.            5.
(1,0,0,*)              3.            3.            3.            3.            3.
(1,0,1,*)              3.            3.            3.            3.            3.
( z6 s' Q0 k; C9 m; n4 r(1,1,0,*)              4.            4.            4.            4.            4.
(1,1,1,*)              4.            4.            4.            4.            4.
(1,2,0,*)              5.            5.            5.            5.            5.
(1,2,1,*)              5.            5.            5.            5.            5.
+ n9 b6 n) F' c+ F4 S4 c+ w! J1 Y
& Z) g: a- }8 I8 k(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.
(0,2,0,*)              6.            6.            6.            6.            6.
(0,2,1,*)              7.            7.            7.            7.            7.
(1,0,0,*)              6.            6.            6.            6.            6.
(1,0,1,*)              7.            7.            7.            7.            7.
. _6 U9 J/ G4 H(1,1,0,*)              6.            6.            6.            6.            6.
(1,1,1,*)              7.            7.            7.            7.            7.
$ [5 R1 E, t8 |2 T" ~9 G(1,2,0,*)              6.            6.            6.            6.            6.
1 o5 |" K* f" z(1,2,1,*)              7.            7.            7.            7.            7.

(0,0,0,*)              8.            9.           10.           11.           12.
0 h2 o7 Q5 K* H' {8 p(0,0,1,*)              8.            9.           10.           11.           12.
: Z: P0 q$ g' ]2 E3 B$ A2 {' N2 _(0,1,0,*)              8.            9.           10.           11.           12.
5 Q( q1 G) C+ U- w(0,1,1,*)              8.            9.           10.           11.           12.
8 v7 h/ _3 l; T( ~( X" j; Y$ i7 G; W(0,2,0,*)              8.            9.           10.           11.           12.
(0,2,1,*)              8.            9.           10.           11.           12.
5 T" n+ Q. M. k- E3 \' m(1,0,0,*)              8.            9.           10.           11.           12.
' M0 y$ W: e6 H(1,0,1,*)              8.            9.           10.           11.           12.
(1,1,0,*)              8.            9.           10.           11.           12.
  d' m, T, t( I3 O% n(1,1,1,*)              8.            9.           10.           11.           12.
(1,2,0,*)              8.            9.           10.           11.           12.
(1,2,1,*)              8.            9.           10.           11.           12.
/ }' [  e' Q9 ^4 z0 c
% K. f. U& @; O8 q/ V* W(0,0,0,*)             18.           19.           20.           21.           22.
8 F. Y- {+ u( a(0,0,1,*)             19.           20.           21.           22.           23.
+ u8 C$ \1 U* D1 M' O' i" A  }; n& i(0,1,0,*)             19.           20.           21.           22.           23.
' w. R1 w" X2 ]6 ]' J, c- S(0,1,1,*)             20.           21.           22.           23.           24.
(0,2,0,*)             20.           21.           22.           23.           24.
8 N7 R$ T2 h. x(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 G7 g4 |4 R0 A: z& [# X(1,1,0,*)             20.           21.           22.           23.           24.
# I3 ^; V7 p1 ]8 N- S(1,1,1,*)             21.           22.           23.           24.           25.
) R# f5 W9 x" B. ?' f; j; |(1,2,0,*)             21.           22.           23.           24.           25.
( Y* D! K( C0 k* O& j8 l(1,2,1,*)             22.           23.           24.           25.           26.

(0,0,*)            100.          105.
(0,1,*)            105.          110.
; R* T; |4 j) z  l3 V  A- m(0,2,*)            110.          115.
0 }- m* Y( z# g$ J1 F(1,0,*)            105.          110.
; E. C$ Q- k! N5 C; M5 x(1,1,*)            110.          115.
2 k' z, _9 f* [  P, H  a(1,2,*)            115.          120.

(0,*)            205.          215.          225.
6 Q0 e$ P( D! U3 r# b(1,*)            215.          225.          235.

(0,*)            645.
1 q7 q, x8 w4 E8 y0 h(1,*)            675.

& w: U/ d. l- q1320.

; I1 t: @6 d& G; d6 r6 e

forcal

在matlab中，纯for循环速度最慢。而一半for循环+一半向量化的速度最快，Forcal中也是如此：
) n) j: T3 x; Q5 D% V7 N* D7 y

2. !using["math","sys"];\\" V+ ^( R8 [% ~2 d. c
   
4. mvar:0 W* B' u8 \; X
   
6. (:p1,p2,p3,a,b)=
   
7. ; H0 F4 O+ t0 W7 G% ^5 z
8. {, d7 }! P8 q& i6 c
   
10.   oo{
    
11. 4 q( O: H( @4 T\\" [3 }1 O
12.     a=array[1000].rand(),. P/ d, g! N. y; e/ R
    
14.     b=array[1000].rand(),
    
15. ) ^, A* G% R4 l
16.     p1=array[1000,1000],
    
17. + C$ r: w3 Q, }
18.     p2=array[1000,1000],
    
19. 0 C2 a7 c' K5 A5 ^; a! j
20.     p3=array[1000,1000],; ?+ h0 R# {. E9 H
    
22.     t0=clock(),
    
23. . W0 E+ O7 y/ r( g- X
24.     ndgrid(a,b,&A,&B),
    
25. . U  S2 \- y6 f\\" `5 F
26.     p1.=A+B
    
27. ) W. H' X9 `- n0 D( ~- M! x9 I7 T
28.   },; t: v* @( Y' h! W/ o; s  o3 f
    
30.   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},
    
31. , H, W6 A0 h! P7 A
32.   lena=FCDLen(a),
    
33. \\" P' [$ {! o6 {. p; s\\" E& @
34.   lenb=FCDLen(b),
    
35. 5 u: B1 H! C9 W* s  ^& r
36.   t0=clock(),
    
37. $ `2 Z3 N  B+ A- Y- [
38.   m = lenb-1, (m>=0).while{
    
39. + h: G\\" K# [6 z* \* J
40.     oo{p2(m,neg) = a+rn[b(m)]},6 o# a# {& `3 N, l0 `% @. f
    
42.     m--
    
43. * H' f* \: C9 @& s& }( o
44.   },8 a9 c3 s8 C6 A\\" V- x/ F! y\\" I
    
46.   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},% w4 I% J: a& i! U1 ^
    
48.   t0=clock(),; z& U2 f5 a& B+ v8 L! e  j
    
50.   m = lenb-1, (m>=0).while{9 o% l& m  A- S1 _; Z% X
    
52.     n = lena-1, (n>=0).while{
    
53. 3 |$ j, N9 U9 T. c8 b% s$ e
54.       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些
    
55. 6 Q3 J- P+ j- e# a5 U* i\\" h5 b
56.       A(p3,m,n) = A(a,n)+A(b,m),
    
57. \\" A: t; k4 E: o# C
58.       n--
    
59. 6 r! l! D9 C2 z8 ~% V; w% k5 u
60.     },: \- Q7 O( G% S\\" V; q3 Y: R5 {
    
62.     m--/ k9 ^- p  U1 f5 L
    
64.   },( @6 `4 m0 [1 E( _' F1 r% M1 M% E3 Y
    
66.   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}+ y$ C, v! f7 q9 {, n3 @\\" }. @3 a
    
68. };4 H; X0 U$ L& _' k) b, S6 G$ p, @7 j
    

!using["math","sys"];<font class=\\"jammer\\">\\" V+ ^( R8 [% ~2 d. c</font><br /> mvar:<font class=\\"jammer\\">0 W* B' u8 \; X</font><br /> (:p1,p2,p3,a,b)=<br /> <span style=\\"display:none\\">; H0 F4 O+ t0 W7 G% ^5 z</span> {<font class=\\"jammer\\">, d7 }! P8 q& i6 c</font><br />   oo{<br /> <span style=\\"display:none\\">4 q( O: H( @4 T\\" [3 }1 O</span>     a=array[1000].rand(),<font class=\\"jammer\\">. P/ d, g! N. y; e/ R</font><br />     b=array[1000].rand(),<br /> <span style=\\"display:none\\">) ^, A* G% R4 l</span>     p1=array[1000,1000],<br /> <span style=\\"display:none\\">+ C$ r: w3 Q, }</span>     p2=array[1000,1000],<br /> <span style=\\"display:none\\">0 C2 a7 c' K5 A5 ^; a! j</span>     p3=array[1000,1000],<font class=\\"jammer\\">; ?+ h0 R# {. E9 H</font><br />     t0=clock(),<br /> <span style=\\"display:none\\">. W0 E+ O7 y/ r( g- X</span>     ndgrid(a,b,&A,&B),<br /> <span style=\\"display:none\\">. U  S2 \- y6 f\\" `5 F</span>     p1.=A+B<br /> <span style=\\"display:none\\">) W. H' X9 `- n0 D( ~- M! x9 I7 T</span>   },<font class=\\"jammer\\">; t: v* @( Y' h! W/ o; s  o3 f</font><br />   printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},<br /> <span style=\\"display:none\\">, H, W6 A0 h! P7 A</span>   lena=FCDLen(a),<br /> <span style=\\"display:none\\">\\" P' [$ {! o6 {. p; s\\" E& @</span>   lenb=FCDLen(b),<br /> <span style=\\"display:none\\">5 u: B1 H! C9 W* s  ^& r</span>   t0=clock(),<br /> <span style=\\"display:none\\">$ `2 Z3 N  B+ A- Y- [</span>   m = lenb-1, (m>=0).while{<br /> <span style=\\"display:none\\">+ h: G\\" K# [6 z* \* J</span>     oo{p2(m,neg) = a+rn[b(m)]},<font class=\\"jammer\\">6 o# a# {& `3 N, l0 `% @. f</font><br />     m--<br /> <span style=\\"display:none\\">* H' f* \: C9 @& s& }( o</span>   },<font class=\\"jammer\\">8 a9 c3 s8 C6 A\\" V- x/ F! y\\" I</font><br />   printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},<font class=\\"jammer\\">% w4 I% J: a& i! U1 ^</font><br />   t0=clock(),<font class=\\"jammer\\">; z& U2 f5 a& B+ v8 L! e  j</font><br />   m = lenb-1, (m>=0).while{<font class=\\"jammer\\">9 o% l& m  A- S1 _; Z% X</font><br />     n = lena-1, (n>=0).while{<br /> <span style=\\"display:none\\">3 |$ j, N9 U9 T. c8 b% s$ e</span>       //p3(m,n) = a(n)+b(m),  //用这句还要慢一些<br /> <span style=\\"display:none\\">6 Q3 J- P+ j- e# a5 U* i\\" h5 b</span>       A(p3,m,n) = A(a,n)+A(b,m),<br /> <span style=\\"display:none\\">\\" A: t; k4 E: o# C</span>       n--<br /> <span style=\\"display:none\\">6 r! l! D9 C2 z8 ~% V; w% k5 u</span>     },<font class=\\"jammer\\">: \- Q7 O( G% S\\" V; q3 Y: R5 {</font><br />     m--<font class=\\"jammer\\">/ k9 ^- p  U1 f5 L</font><br />   },<font class=\\"jammer\\">( @6 `4 m0 [1 E( _' F1 r% M1 M% E3 Y</font><br />   printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}<font class=\\"jammer\\">+ y$ C, v! f7 q9 {, n3 @\\" }. @3 a</font><br /> };<font class=\\"jammer\\">4 H; X0 U$ L& _' k) b, S6 G$ p, @7 j</font><br />

结果：
ndgrid: 3.2001e-002
for1: 1.4999e-002
for2: 1.86

1 L! U1 ^1 s2 u( w& o

forcal

一段程序的Forcal实现：
3 T& V( H% d: y* Y2 P2 N+ c! X
: ^* V3 q. N. D. r//用C++代码描述为：
s=0.0;
for(x=0.0;x<=1.0;x=x+0.0011)
1 j! W% a0 |% {$ R5 b7 k1 M{
0 ~+ @# J) F1 S, Y0 ?/ Z# C$ _* O   for(y=1.0;y<=2.0;y=y+0.0009)
* H* F. d; b% m& Z1 A   {
7 A7 q( [2 g7 c0 q4 @9 a- h0 P/ p     s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   }
" @* e- Q1 t5 n" K  |* `& R: m}

1、**数组求和函数Sum

* T: R& s4 ]0 b

2. !using["math","sys"];/ ~1 \6 C7 ?7 r, L# g
   
4. mvar:
   
5. ! N: L- t2 B% [! T3 a/ D
6. t=clock(),5 i+ y8 n6 Q2 w; r/ e\\" v  I4 _
   
8. oo{0 v\\" b, y+ e: s/ Y  c4 Y9 Y7 G$ m
   
10.   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
    
11. & H9 z4 O; i- C5 l/ F- S1 w
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. 0 v! i, u7 c\\" p7 @3 i2 A5 L, l
14. };7 {/ Y2 d# h5 Z3 ~' R
    
16. [clock()-t]/1000;

!using["math","sys"];<font class=\\"jammer\\">/ ~1 \6 C7 ?7 r, L# g</font><br /> mvar:<br /> <span style=\\"display:none\\">! N: L- t2 B% [! T3 a/ D</span> t=clock(),<font class=\\"jammer\\">5 i+ y8 n6 Q2 w; r/ e\\" v  I4 _</font><br /> oo{<font class=\\"jammer\\">0 v\\" b, y+ e: s/ Y  c4 Y9 Y7 G$ m</font><br />   ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],<br /> <span style=\\"display:none\\">& H9 z4 O; i- C5 l/ F- S1 w</span>   Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn  (2)))))),0]<br /> <span style=\\"display:none\\">0 v! i, u7 c\\" p7 @3 i2 A5 L, l</span> };<font class=\\"jammer\\">7 {/ Y2 d# h5 Z3 ~' R</font><br /> [clock()-t]/1000;

! X( J4 Q5 v# t1 |; ~5 o结果：
1008606.64947441
. ^" ~, W9 \! J) N3 ]( I0.625   //时间

. P8 m4 I" s8 @& U2、求和函数sum
  }4 \- ^  p2 q; C/ a8 ~
/ v3 u0 V) D5 @  z

1. f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2)))));
   % F- s# w. G+ q\" f! P
2. sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
( [5 y3 k  Z/ h& E" U# m- w1008606.64947441
0.719   //时间

2 Q/ r, S8 Z5 w3、while循环
& l! R  ^- t5 p  A/ H5 K& v7 Y: ~
0 m  B; M4 U; Z8 G

1. mvar:\" ]) \3 J+ n9 Z( {9 Z
   
2. t=sys::clock();( O4 h# ?3 }8 X7 X. E5 ]& V5 ~\" k. `# R
   
3. s=0,x=0,
   . ?2 G5 h, v9 r7 h9 ~$ Q
4. while{x<=1,  //while循环算法；   {. N8 P7 s0 u4 [( N
   
5.    y=1,
   : m  F6 h* X4 x1 b
6.    while{y<=2,
   + G  Y2 p* T, G
7.        s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
   4 ~# N& p! l( i# ?. m! [$ K
8.        y=y+0.0009
   # n- G( `: B1 [0 J: C/ C
9.       },
   \" E% f# v\" s5 }' j; t5 ~! O
10.    x=x+0.0011 1 d1 [4 \' M* _7 a2 @' ?
    
11. }, $ \9 C: x' i, W6 x' M5 `7 w$ m
    
12. s;# O; C% N! f. H$ p0 H
    
13. [sys::clock()-t]/1000;

复制代码

结果：
8 S1 P8 v# r7 o2 y1 }* G1008606.64947441
0.734   //时间

qbist

好深奥！~~~~

forcal

. V/ r; X" P/ V. ^) r0 d) `

好深奥！~~~~
qbist 发表于 2010-10-7 14:56

7 ^2 p' r% e9 W- ]先了解一下，以备不时之需，有问题可以交流哦，呵呵。

qbist

回复 forcal 的帖子
* }7 d" Z% b0 S8 A. X8 ?4 `& y

# R0 d/ v4 k1 K  m. P) E, T. ~3 a' p. E. b    嗯！！！

forcal

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

以下是FcMath中的矩阵乘与徐士良算法库XSLSF（普通的C/C++算法）中的矩阵乘的效率比较：

3 w3 [  a& C8 V: ^1、FcMath中的矩阵乘

2. !using["math","sys"];0 }' y4 G9 k: a, g; m! m7 a* Z& x
   
4. (:a,b,k,t0)=7 d/ @8 Q, k' J0 U8 G9 k9 S
   
6. oo{
   
7. 4 g& Q% i& c$ k9 m; O. M7 o; S
8.   a=rand[1000,1000], b=rand[1000,1000],
   
9. / }4 Z: J, t. C
10.   t0=clock(),2 }- r& O1 u& A8 n8 e  _4 L
    
12.   k=a*b,  //矩阵乘
    
13. + l8 J' ^9 }8 d
14.   k[1,3:5,9].outm()8 i8 S) g+ I) i% ^& O
    
16. },) m8 {' d+ U, h- ^) P9 }7 s. i- u2 t
    
18. [clock()-t0]/1000;* F$ a2 P. V4 V* K+ O! W$ d/ x
    

!using["math","sys"];<font class=\\"jammer\\">0 }' y4 G9 k: a, g; m! m7 a* Z& x</font><br /> (:a,b,k,t0)=<font class=\\"jammer\\">7 d/ @8 Q, k' J0 U8 G9 k9 S</font><br /> oo{<br /> <span style=\\"display:none\\">4 g& Q% i& c$ k9 m; O. M7 o; S</span>   a=rand[1000,1000], b=rand[1000,1000],<br /> <span style=\\"display:none\\">/ }4 Z: J, t. C</span>   t0=clock(),<font class=\\"jammer\\">2 }- r& O1 u& A8 n8 e  _4 L</font><br />   k=a*b,  //矩阵乘<br /> <span style=\\"display:none\\">+ l8 J' ^9 }8 d</span>   k[1,3:5,9].outm()<font class=\\"jammer\\">8 i8 S) g+ I) i% ^& O</font><br /> },<font class=\\"jammer\\">) m8 {' d+ U, h- ^) P9 }7 s. i- u2 t</font><br /> [clock()-t0]/1000;<font class=\\"jammer\\">* F$ a2 P. V4 V* K+ O! W$ d/ x</font><br />

结果：

1.         238.447        247.837        247.065        248.105        247.0583 a; ^# V3 S# I\" Z
   
2.         244.123        249.925        247.553        243.981        250.016
   ' Z8 K\" K0 z\" z$ X% P& t/ `
3.         236.387        252.025        245.651        248.866        248.866
     d0 O& l9 v+ V, w& C) [
4. 2.219 秒/ y% x& a$ S# {9 ^. _. o3 `5 E$ G
   

复制代码

) S( I, \9 W. |& ?% V  z; b5 ?* K' m, t2、XSLSF（普通的C/C++算法）中的矩阵乘

2. !using["math","sys","XSLSF"];
   
3. 4 m9 M9 `$ O% A- t
4. (:a,b,k,t0)=( o4 D9 O- |$ x\\" {, ~) P
   
6. oo{
   
7. 9 o6 C9 e) C) [* }% W
8.   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],, B5 z9 D3 f9 A+ @& Q, \
   
10.   t0=clock(),
    
11. 6 A2 G- J) _  Y3 H: {* ?$ |
12.   rmul[k:a,b],  //矩阵乘, {) P/ C8 Z) j. s, ^3 _$ @
    
14.   k[1,3:5,9].outm()
    
15. 6 w; |: V; ^) H
16. },
    
17. 3 |1 @% W3 e% z\\" p* E: R# o
18. [clock()-t0]/1000;
    
19. + [5 n' E6 x5 r7 Q

!using["math","sys","XSLSF"];<br /> <span style=\\"display:none\\">4 m9 M9 `$ O% A- t</span> (:a,b,k,t0)=<font class=\\"jammer\\">( o4 D9 O- |$ x\\" {, ~) P</font><br /> oo{<br /> <span style=\\"display:none\\">9 o6 C9 e) C) [* }% W</span>   a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],<font class=\\"jammer\\">, B5 z9 D3 f9 A+ @& Q, \</font><br />   t0=clock(),<br /> <span style=\\"display:none\\">6 A2 G- J) _  Y3 H: {* ?$ |</span>   rmul[k:a,b],  //矩阵乘<font class=\\"jammer\\">, {) P/ C8 Z) j. s, ^3 _$ @</font><br />   k[1,3:5,9].outm()<br /> <span style=\\"display:none\\">6 w; |: V; ^) H</span> },<br /> <span style=\\"display:none\\">3 |1 @% W3 e% z\\" p* E: R# o</span> [clock()-t0]/1000;<br /> <span style=\\"display:none\\">+ [5 n' E6 x5 r7 Q</span>

结果：

1.         262.121        247.583        260.529        259.548        258.328% u4 r8 F& [. a
   
2.         255.413        246.563        254.356        250.548        251.509' z% d* u+ G0 H7 P5 B6 c3 @. F
   
3.         256.152        247.725        259.444        250.827        249.8168 l5 D! u6 g% D, }1 ^; D) v$ r* |1 f6 L
   
4. 10.563 秒& _8 z% H2 a* }0 s# N+ H
   

复制代码