查看: 5727|回复: 9

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

[复制链接]

字体大小: 正常放大

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

电梯直达

1^#

发表于 2010-10-7 11:45 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

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！如果您有什么好的算法，任何改进的意见或建议，请与作者联系。

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

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

zan

已有 1 人评分	体力	收起理由
厚积薄发	+ 5

总评分: 体力 + 5 查看全部评分

转播0 淘帖0 分享0 收藏0 支持0 反对0 微信

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

2^#

发表于 2010-10-7 11:48 |只看该作者 |招呼Ta 关注Ta

例子1代码：

!using["math"];, J# @5 j: Y* ?
mvar:
1 E- ]! @; R% c; e- x3 I. k
oo{ //一般在oo函数中调用FcMath函数% Q! a8 U1 ^0 _ g: K5 ^
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化
, |4 Q( z! K. r
a.outm(), //输出矩阵a u\\" _: d; i) v\\" v
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出! Y+ s) P8 e. x( F
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出+ w) m& M( [3 E& C j& k0 N5 s1 F
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出
3 t+ I8 U- L- z2 N' t
};: V# y4 ?( h6 P/ G7 N! I3 }\\" G

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.852615+ ]+ G- h! ], ` K6 s
0.630646 0.927048 0.440308 0.162857 0.556854
% q% W$ o* ` Z- u
0.43309 0.34552 0.563919 0.937164 0.209641
5 U9 z& x/ T$ Y\" K, T9 @, h) ?
0.603271 0.727676 0.130951 5.35736e-002 0.197937
3 U2 o* h. }' J5 M
0.576004 0.747589 1.17645e-002 0.363892 0.280777
$ j5 ^. f0 D- D( P
0.646454 0.381088 0.58551 0.26387 0.93692
+ k# ^& N* F# A0 H
4 D0 {* Y; ~; L; B, a) e* x
0.153259' m. p- p, E& S9 n
0.162857; @+ k3 j i! R4 Y7 L) U
0.937164* [- ^( Q, X1 X3 |# c! J
5.35736e-002: U# p% o2 J& f/ V2 O
0.363892: ]. f2 @( O\" e
0.263871 B+ h6 S S* h. J4 B3 I
# \6 T: W( \- R* e: S$ u
0.603271 0.727676 0.130951 5.35736e-002 0.197937
; j8 ~: ~( {2 a4 ?
* g5 U/ a+ o( A7 y- q5 S9 k
0.130951 5.35736e-002
$ X+ ?5 Y3 ?% ~4 \! _6 K3 ]% [1 e# V$ X
1.17645e-002 0.3638923 t- X$ H\" T+ }
0.58551 0.263872 V3 y! u, Y8 q; X3 t8 P' a
O+ B% T( S% @

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义
A; Y: s0 c6 |
{' R8 e [, \% {, M2 u- Z& w
y1=x1*x1+x2*x2+x3*x3-1.0,: z! c$ Z9 e, N
y2=2.0*x1*x1+x2*x2-4.0*x3,' F\\" @0 n# d1 Y) ~8 k# Q
y3=3.0*x1*x1-4.0*x2+x3*x3! a8 E' T& h) X# @+ r
}; g% a& I* \, X% L
!using["math","sys"];
) |5 s% k9 a, D5 N1 `
mvar:7 T' k) y3 A/ O: s7 x5 ~
oo{5 h, V4 T0 ?: W2 q
x=array(3),
; }. Z6 y( m: K8 ?* v6 ]+ s- S
x.SA[0 : 1,1,1], //设置初值为1,1,11 w' g5 X: Y$ e
i=netn[HFor("f"),x], //拟牛顿法解方程! C3 }9 s% p) |0 ]& J
x.outm(), //输出结果
( d6 f' D' N$ n
i //返回迭代次数& y& F3 f* O/ F, R# C8 g
};' g+ a+ }5 C# f- b

f(x1,x2,x3,y1,y2,y3)=      //函数定义 
  A; Y: s0 c6 |
{' R8 e  [, \% {, M2 u- Z& w

y1=x1*x1+x2*x2+x3*x3-1.0,: z! c$ Z9 e, N

y2=2.0*x1*x1+x2*x2-4.0*x3,' F\\" @0 n# d1 Y) ~8 k# Q

y3=3.0*x1*x1-4.0*x2+x3*x3! a8 E' T& h) X# @+ r

};  g% a& I* \, X% L

!using["math","sys"]; 
) |5 s% k9 a, D5 N1 `
mvar:7 T' k) y3 A/ O: s7 x5 ~

oo{5 h, V4 T0 ?: W2 q

x=array(3), 
; }. Z6 y( m: K8 ?* v6 ]+ s- S
  x.SA[0 : 1,1,1],       //设置初值为1,1,11 w' g5 X: Y$ e

i=netn[HFor("f"),x],   //拟牛顿法解方程! C3 }9 s% p) |0 ]& J

x.outm(),              //输出结果 
( d6 f' D' N$ n
  i                      //返回迭代次数& y& F3 f* O/ F, R# C8 g

};' g+ a+ }5 C# f- b

结果：
0.785197 0.496611 0.369923

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

3^#

发表于 2010-10-7 11:53 |只看该作者 |招呼Ta 关注Ta

效率测试：
simwe的网友lin2009 的matlab代码：

clear all
) w+ l7 J! w9 q
clc
k) n1 b) z1 ^) i* @
tic
9 M! C4 Q9 n/ b, D/ ^
k = zeros(5,5); % //生成5×5全0矩阵
, E1 {9 S- @0 s$ f
% 循环计算以下程序段100000次：
4 {4 w5 ^8 L% c. N C
for m = 1:100000$ X' r, C; x: v% k\" M6 @
a = rand(5,7);
/ Z. y# P% D\" O( T6 P) V
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化5 L* Z- e. q. I1 m% z4 }3 r
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);( k1 m d, {7 ~- p# T7 ~\" M
end
4 y3 q8 v) G8 L1 j4 v
k
) s/ E$ U# R' I& d1 r; a
toc\" y6 y7 ^7 E; @0 u9 ?7 G

复制代码

Forcal代码：

运行稍快的代码，比matlab约快10％吧？

!using["math","sys"];$ O! r1 l: t, ]! f, F
mvar:
2 T% U- K4 U, ]6 u) b
t0=clock(),
+ `% @$ Q) t* U1 ]: c; l( I. @
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0* M. ]* q( W2 c1 k- {8 U; T
i=0,(i<1000 00).while{ //循环计算1000 00次
1 R( O- p) b: U& \0 n3 D+ a% ?9 m
oo{
' E! z# u5 V/ V
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化 A) b% z, Z\\" U, O\\" s\\" o
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) [. p5 B/ A' n, L# ~/ ?
},
5 a\\" i$ t# x; ]! y
i++$ L! ^1 s b U% ^7 ^: d
},0 n0 G8 W' ^+ E
k.outm(), //输出矩阵k，然后销毁k
4 F8 i: S( M, G# K
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"];$ O! r1 l: t, ]! f, F

mvar: 
2 T% U- K4 U, ]6 u) b
t0=clock(), 
+ `% @$ Q) t* U1 ]: c; l( I. @
oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0* M. ]* q( W2 c1 k- {8 U; T

i=0,(i<1000 00).while{ //循环计算1000 00次 
1 R( O- p) b: U& \0 n3 D+ a% ?9 m
  oo{ 
' E! z# u5 V/ V
    a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化  A) b% z, Z\\" U, O\\" s\\" o

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)  [. p5 B/ A' n, L# ~/ ?

}, 
5 a\\" i$ t# x; ]! y
  i++$ L! ^1 s  b  U% ^7 ^: d

},0 n0 G8 W' ^+ E

k.outm(),             //输出矩阵k，然后销毁k 
4 F8 i: S( M, G# K
[clock()-t0]/1000;    //得到计算时间，秒

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

比较好看些的代码，似乎也比matlab稍快吧？

!using["math","sys"];# X# L. m: R* u) x/ v
(:t0,k,i,a,b)=! K' u, e1 o. Y1 d
{/ M\\" d\\" j8 B4 J x' D% D3 L& l
t0=clock(),
+ v! J. V% g5 p# [* V
k=zeros[5,5],1 n2 J0 t% N8 _0 x$ L
i=0,(i<1000 00).while{) V$ D5 Q- |& D$ @, j, |+ Z
oo{/ i5 s4 B- h* S# K/ a# G
a=rand[5,7], b=rand[7,5],; Q3 b# ?\\" t) G# k' m0 \8 K/ ?
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)8 F% p( I1 \) K$ e3 a7 ~6 l
},, ~! H2 Z/ B7 I. Z. o
i++ ]+ i N7 o; B# {/ S- u
},# \5 n$ b3 g+ ^$ |: ], n: W$ {
k.outm().delete(),
& o0 V: x. e) a1 ?% ^
[clock()-t0]/1000) t# y% Z# u5 o0 n6 `
};

!using["math","sys"];# X# L. m: R* u) x/ v

(:t0,k,i,a,b)=! K' u, e1 o. Y1 d

{/ M\\" d\\" j8 B4 J  x' D% D3 L& l

t0=clock(), 
+ v! J. V% g5 p# [* V
  k=zeros[5,5],1 n2 J0 t% N8 _0 x$ L

i=0,(i<1000 00).while{) V$ D5 Q- |& D$ @, j, |+ Z

oo{/ i5 s4 B- h* S# K/ a# G

a=rand[5,7], b=rand[7,5],; Q3 b# ?\\" t) G# k' m0 \8 K/ ?

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)8 F% p( I1 \) K$ e3 a7 ~6 l

},, ~! H2 Z/ B7 I. Z. o

i++  ]+ i  N7 o; B# {/ S- u

},# \5 n$ b3 g+ ^$ |: ], n: W$ {

k.outm().delete(), 
& o0 V: x. e) a1 ?% ^
  [clock()-t0]/1000) t# y% Z# u5 o0 n6 `

};

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

该例子的理论结果是每个元素均为275000。

我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

4^#

发表于 2010-10-7 11:55 |只看该作者 |招呼Ta 关注Ta

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

!using["math"];
/ H, j( W% [+ x
mvar:
2 U9 }) k$ }6 Z( z8 D
oo{/ E$ ]) w8 ]6 ^, F* H& r
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
, \3 K1 a3 f$ h: b0 y
a1.outm[5,1,1],
2 B2 q4 d- h$ C
a2.outm[5,1,1],: x) _9 h& U7 _1 L8 q- o3 [
a3.outm[5,1,1],1 x, Z& n0 V2 [5 i y
a4.outm[5,1,1],
% l. h! C, S% E; p\\" I
a=a1+a2+a3+a4,
9 \. ~+ k6 P0 m r: l\\" [ K6 S. |
a.outm[5,1,1],: @4 j$ q, z3 t4 v7 ^# A! t7 x
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]2 _5 r: ~& o$ `# z4 e( g1 N7 F% c1 S5 L
};# z- Q3 K! h: B- [. e

!using["math"]; 
/ H, j( W% [+ x
mvar: 
2 U9 }) k$ }6 Z( z8 D
oo{/ E$ ]) w8 ]6 ^, F* H& r

ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4], 
, \3 K1 a3 f$ h: b0 y
  a1.outm[5,1,1], 
2 B2 q4 d- h$ C
  a2.outm[5,1,1],: x) _9 h& U7 _1 L8 q- o3 [

a3.outm[5,1,1],1 x, Z& n0 V2 [5 i  y

a4.outm[5,1,1], 
% l. h! C, S% E; p\\" I
  a=a1+a2+a3+a4, 
9 \. ~+ k6 P0 m  r: l\\" [  K6 S. |
  a.outm[5,1,1],: @4 j$ q, z3 t4 v7 ^# A! t7 x

Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]2 _5 r: ~& o$ `# z4 e( g1 N7 F% c1 S5 L

};# z- Q3 K! h: B- [. e

说明：
linspace(8,12,5)：生成一维数组，共5个元素8～12
a1.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）：

(0,0,0,*) 1. 1. 1. 1. 1.
(0,0,1,*) 1. 1. 1. 1. 1.
(0,1,0,*) 1. 1. 1. 1. 1.
(0,1,1,*) 1. 1. 1. 1. 1.
(0,2,0,*) 1. 1. 1. 1. 1.
(0,2,1,*) 1. 1. 1. 1. 1.
(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.

(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.
(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.

(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.
(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.
(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.
(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.
(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.
(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.
(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.
(1,2,0,*) 21. 22. 23. 24. 25.
(1,2,1,*) 22. 23. 24. 25. 26.

(0,0,*) 100. 105.
(0,1,*) 105. 110.
(0,2,*) 110. 115.
(1,0,*) 105. 110.
(1,1,*) 110. 115.
(1,2,*) 115. 120.

(0,*) 205. 215. 225.
(1,*) 215. 225. 235.

(0,*) 645.
(1,*) 675.

1320.

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

5^#

发表于 2010-10-7 11:56 |只看该作者 |招呼Ta 关注Ta

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

!using["math","sys"];5 v+ Q d1 l- `2 t
mvar:
! z' O) p9 b/ Y- l
(:p1,p2,p3,a,b)=
Y% P\\" P5 s7 R4 s. ^) n
{
- P/ w% m) l! u3 R\\" G8 c
oo{
+ L l/ s9 r. i0 p( J
a=array[1000].rand(),
7 n, ~1 C8 w% F* L# T
b=array[1000].rand(),1 f; A, @: L0 x3 E3 g/ L9 }
p1=array[1000,1000],
9 m' j$ ^# R Z! v* U3 c
p2=array[1000,1000],/ G& |\\" g3 u# G5 u& W
p3=array[1000,1000],5 D4 L! w( t4 l/ `, b2 [8 t
t0=clock(),2 `+ p5 o6 ]4 z+ V9 P8 ]
ndgrid(a,b,&A,&B),0 ^4 t7 V( B( ?
p1.=A+B$ M2 D3 E* `' o' r4 r, C
},
0 r0 u+ c u' N) y5 f' k( Y: F+ P+ }
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},
1 J- R2 f& F8 \1 _6 W( I
lena=FCDLen(a),& C& k. a/ n! j D; i
lenb=FCDLen(b),' G9 o) G' C( m: t# l# j
t0=clock(),# s6 c1 K$ k4 I# z: X
m = lenb-1, (m>=0).while{( @ D+ M* U8 C\\" M& |
oo{p2(m,neg) = a+rn[b(m)]},
\\" V5 u. ? R' g+ a- F, o5 C) a
m--( C2 v2 h% I, j ^2 l1 e
},8 `4 G1 w- I- \5 \& C2 n/ V( K
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},+ a& c; ~! V+ g5 V# |' J
t0=clock(),1 H$ Y; Z2 }9 S, t
m = lenb-1, (m>=0).while{
' c) }7 b' U\\" j# `\\" o
n = lena-1, (n>=0).while{
! G( q8 W T\\" |4 }0 X+ Z9 c4 M
//p3(m,n) = a(n)+b(m), //用这句还要慢一些
0 a; `) I5 i% H+ x
A(p3,m,n) = A(a,n)+A(b,m),
3 c' W6 C& K/ C$ `. ?
n--
! a {, b7 d4 a; b5 p; ^& ~ u- v
},$ ^0 u& _9 }% {; N' j9 f
m--6 o% k6 z) O! @# \
},3 d( i% \5 T% M$ R
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}3 X3 o0 B9 E3 P4 N& ]2 H
};4 _; ^, P. ?6 ^. a! ]* |

!using["math","sys"];5 v+ Q  d1 l- `2 t

mvar: 
! z' O) p9 b/ Y- l
(:p1,p2,p3,a,b)= 
  Y% P\\" P5 s7 R4 s. ^) n
{ 
- P/ w% m) l! u3 R\\" G8 c
  oo{ 
+ L  l/ s9 r. i0 p( J
    a=array[1000].rand(), 
7 n, ~1 C8 w% F* L# T
    b=array[1000].rand(),1 f; A, @: L0 x3 E3 g/ L9 }

p1=array[1000,1000], 
9 m' j$ ^# R  Z! v* U3 c
    p2=array[1000,1000],/ G& |\\" g3 u# G5 u& W

p3=array[1000,1000],5 D4 L! w( t4 l/ `, b2 [8 t

t0=clock(),2 `+ p5 o6 ]4 z+ V9 P8 ]

ndgrid(a,b,&A,&B),0 ^4 t7 V( B( ?

p1.=A+B$ M2 D3 E* `' o' r4 r, C

}, 
0 r0 u+ c  u' N) y5 f' k( Y: F+ P+ }
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000}, 
1 J- R2 f& F8 \1 _6 W( I
  lena=FCDLen(a),& C& k. a/ n! j  D; i

lenb=FCDLen(b),' G9 o) G' C( m: t# l# j

t0=clock(),# s6 c1 K$ k4 I# z: X

m = lenb-1, (m>=0).while{( @  D+ M* U8 C\\" M& |

oo{p2(m,neg) = a+rn[b(m)]}, 
\\" V5 u. ?  R' g+ a- F, o5 C) a
    m--( C2 v2 h% I, j  ^2 l1 e

},8 `4 G1 w- I- \5 \& C2 n/ V( K

printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},+ a& c; ~! V+ g5 V# |' J

t0=clock(),1 H$ Y; Z2 }9 S, t

m = lenb-1, (m>=0).while{ 
' c) }7 b' U\\" j# `\\" o
    n = lena-1, (n>=0).while{ 
! G( q8 W  T\\" |4 }0 X+ Z9 c4 M
      //p3(m,n) = a(n)+b(m),  //用这句还要慢一些 
0 a; `) I5 i% H+ x
      A(p3,m,n) = A(a,n)+A(b,m), 
3 c' W6 C& K/ C$ `. ?
      n-- 
! a  {, b7 d4 a; b5 p; ^& ~  u- v
    },$ ^0 u& _9 }% {; N' j9 f

m--6 o% k6 z) O! @# \

},3 d( i% \5 T% M$ R

printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}3 X3 o0 B9 E3 P4 N& ]2 H

};4 _; ^, P. ?6 ^. a! ]* |

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

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

6^#

发表于 2010-10-7 11:58 |只看该作者 |招呼Ta 关注Ta

一段程序的Forcal实现：

//用C++代码描述为：
s=0.0;
for(x=0.0;x<=1.0;x=x+0.0011)
{
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)))));
}
}

1、**数组求和函数Sum

!using["math","sys"];: A, G0 U0 \( @' r
mvar:
/ O& u! M6 Y$ S: O; X! t+ b4 u: P
t=clock(),
$ c4 U% S+ G! H0 c. h! o- Z
oo{. n e- n4 y4 }/ f
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],3 u- K6 P8 }$ [/ c* O. O+ V' f/ E
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]- A5 @! T8 V# b% f* ~! M f& D) v
};
, ]( @( K: a8 n* n
[clock()-t]/1000;

结果：
1008606.64947441
0.625 //时间

2、求和函数sum

f(x,y)=cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))); 6 Q+ M2 }+ s3 b3 k( A i5 x; c& O
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:
( } |1 ~' M X. F' |! B% d6 Z
t=sys::clock();
`& T9 k2 s; Q
s=0,x=0, 6 ~1 r! @; ?, ?1 ] `
while{x<=1, //while循环算法；
$ o3 A2 ~/ M5 o, p, k
y=1,
U* w( q- u4 [4 D- w4 R
while{y<=2,
5 y' ^% h1 a4 c8 _; L2 b# J1 w
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 2 Q- w0 h7 o2 J. |# p
y=y+0.0009
0 O7 O1 N. { }& s3 i
},
: G5 L) {# \6 S: ]3 b6 ~* ?
x=x+0.0011 ; i, h* ]) n. C: J# d8 q
},
/ J! U: K2 a# y. B7 k$ x3 r
s;
# C0 g* O c* L- {9 o% O9 f$ ~
[sys::clock()-t]/1000;

复制代码

结果：
1008606.64947441
0.734 //时间

我也说一句

发表

使用道具举报

qbist

2 主题	3 听众	304 积分

升级 1.33%

该用户从未签到

自我介绍: 一个对未来充满信心的阳光型男孩！

7^#

发表于 2010-10-7 14:56 |只看该作者 |招呼Ta 关注Ta

好深奥！~~~~

我也说一句

发表

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

8^#

发表于 2010-10-8 21:09 |只看该作者 |招呼Ta 关注Ta

本帖最后由 forcal 于 2010-10-8 21:10 编辑

好深奥！~~~~4 I) P/ [6 e1 {# i- N2 v/ h. M4 @
qbist 发表于 2010-10-7 14:56

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

我也说一句

发表

使用道具举报

qbist

2 主题	3 听众	304 积分

升级 1.33%

该用户从未签到

自我介绍: 一个对未来充满信心的阳光型男孩！

9^#

发表于 2010-10-9 15:54 |只看该作者 |招呼Ta 关注Ta

回复 forcal 的帖子

嗯！！！

我也说一句

发表

好好学习天天向上！！！

使用道具举报

forcal

45 主题	3 听众	282 积分

升级 91%

TA的每日心情

	难过 2012-8-27 18:22

签到天数: 1 天

[LV.1]初来乍到

10^#

发表于 2010-10-13 18:52 |只看该作者 |招呼Ta 关注Ta

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

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

1、FcMath中的矩阵乘

!using["math","sys"];
\\" y! {& ~: t* @' `( n$ w
(:a,b,k,t0)=
# p+ t& X0 h# [3 X* d- u
oo{/ ~; b7 Y/ U2 M+ J6 h
a=rand[1000,1000], b=rand[1000,1000],
$ f }7 z5 Y- b
t0=clock(),
2 m3 @: o! w2 q* d' _1 G
k=a*b, //矩阵乘% @) k; d- Y* \' r! W0 d4 K
k[1,3:5,9].outm()7 T& M0 L0 @: v. b
},
# r# L' {: d8 H) ?& C5 O9 `$ Y# J
[clock()-t0]/1000;
% a2 n4 |9 r5 K! D3 U0 H1 r5 V

结果：

238.447 247.837 247.065 248.105 247.058
! K: C# {\" c3 ~3 K% D
244.123 249.925 247.553 243.981 250.016
9 Q; @ P5 v% o1 Y! D3 R
236.387 252.025 245.651 248.866 248.866\" e6 [2 H( ^8 G6 N$ G0 I
2.219 秒' A. {' ^) C! Y' Z$ F

复制代码

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

!using["math","sys","XSLSF"];5 A* Q) b' E8 @ A3 M$ N
(:a,b,k,t0)=; C$ b- u& `7 ~& O0 _( X; r
oo{1 m) h; ~6 `- ]! x# Z4 K' d
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],4 P9 P/ a- U& p8 i0 \/ `+ W
t0=clock(),
: i\\" J/ }: W* y7 A. t1 M4 m
rmul[k:a,b], //矩阵乘
\\" A0 H; R6 B4 Z1 C' p: T. D
k[1,3:5,9].outm()3 Z+ e' S& f# ^7 J2 s
},
1 Z( b\\" W, o* W* X5 L
[clock()-t0]/1000;! K: `4 W) D x) W. @* I/ ?

结果：

262.121 247.583 260.529 259.548 258.328( b, m4 P- E+ v+ f7 \* J3 }
255.413 246.563 254.356 250.548 251.509
2 y8 z\" I3 d0 ~) f) J, T5 s
256.152 247.725 259.444 250.827 249.8161 P! O9 @+ a* ?' @% {1 w
10.563 秒1 t6 O) q9 L0 z6 K% N4 F/ e

复制代码

我也说一句

发表

使用道具举报

返回列表

帐号		密码		只需要一步，快速开始		注册地址	找回密码

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

浏览过的版块

新人进步奖

QQ

电话咨询

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务|