查看: 5457|回复: 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"];
% c+ ?! t& P7 T+ R% x2 m* C% O/ H/ d
mvar:
5 f7 T# l; ~. D# o1 Z6 n* R& }( y
oo{ //一般在oo函数中调用FcMath函数9 B D, {+ L! {6 E. l
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化
* C: N8 D7 F* w2 j
a.outm(), //输出矩阵a# L0 K4 r* z8 i: G& Z, @2 r
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出; _# V' D; q2 Q
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出
- X- \8 u: J3 U- \' D: y* I
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出1 e- o2 C- w9 |\\" o' V
};5 F* V3 n% B& M/ W. `

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.852615- e) V( m, u9 f
0.630646 0.927048 0.440308 0.162857 0.556854
4 \( h9 O* y! I2 d
0.43309 0.34552 0.563919 0.937164 0.209641
, m( g8 B5 Q2 w# k, p2 s/ V, b
0.603271 0.727676 0.130951 5.35736e-002 0.1979371 V: T* C) z9 f; O
0.576004 0.747589 1.17645e-002 0.363892 0.280777
3 D5 |\" c4 A6 R0 a: Z4 ]0 C
0.646454 0.381088 0.58551 0.26387 0.93692
/ o2 U2 B& z8 ?5 ]! x- Z4 i0 e
8 j8 D: o# V! o
0.153259
2 S3 {& _0 c; o3 e
0.162857- z; Z h# [- i, r7 A\" K7 C# z
0.9371644 E* `8 E3 ~& D* }) q0 A1 T
5.35736e-002 l$ [6 A O\" L g
0.3638928 @$ M7 j; Q+ g, V$ G
0.26387. |: h6 Z4 `# Q9 A A
0 a' w# U/ h. |7 V, C
0.603271 0.727676 0.130951 5.35736e-002 0.197937; Z2 n0 N) M* \+ v) P) ^
' X) y\" u. U* D' Q\" P: Y
0.130951 5.35736e-002- j7 a) o w- r
1.17645e-002 0.3638920 w, b, F- c$ e& u) `\" N
0.58551 0.26387+ b4 J/ ]: y# w' W/ \' G
2 A/ ~+ S* u, ^1 p k

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义
+ r- v0 C# _: W, x }) j9 M3 X
{* ?5 A7 `6 a' S* J/ k
y1=x1*x1+x2*x2+x3*x3-1.0,
9 L* f! X0 {* ~ @- P) D9 G ~
y2=2.0*x1*x1+x2*x2-4.0*x3,! u/ l2 d j* J8 e& n+ s
y3=3.0*x1*x1-4.0*x2+x3*x3
& @\\" A8 b3 B\\" ?* O\\" `, e
};
. h\\" ~2 R4 J+ ?+ a6 [
!using["math","sys"];
! w7 ]& l: L% d/ l) m2 b. w
mvar:+ o1 g0 ^0 w/ z
oo{
' n) @2 E6 t, c# g( N
x=array(3),
$ M$ P4 `5 H! A
x.SA[0 : 1,1,1], //设置初值为1,1,1
K# W$ B2 \! ^+ Z$ g
i=netn[HFor("f"),x], //拟牛顿法解方程
. w, L+ \+ j0 _0 ^\\" E6 ?! \/ ]! d
x.outm(), //输出结果0 f+ D+ I8 E2 c8 q7 O. w. m
i //返回迭代次数/ S* {) H# Z8 i$ K* ^0 I8 P
};
2 }/ j9 d) {4 j, B' s2 {

f(x1,x2,x3,y1,y2,y3)=      //函数定义 
+ r- v0 C# _: W, x  }) j9 M3 X
{* ?5 A7 `6 a' S* J/ k

y1=x1*x1+x2*x2+x3*x3-1.0, 
9 L* f! X0 {* ~  @- P) D9 G  ~
    y2=2.0*x1*x1+x2*x2-4.0*x3,! u/ l2 d  j* J8 e& n+ s

y3=3.0*x1*x1-4.0*x2+x3*x3 
& @\\" A8 b3 B\\" ?* O\\" `, e
}; 
. h\\" ~2 R4 J+ ?+ a6 [
!using["math","sys"]; 
! w7 ]& l: L% d/ l) m2 b. w
mvar:+ o1 g0 ^0 w/ z

oo{ 
' n) @2 E6 t, c# g( N
  x=array(3), 
$ M$ P4 `5 H! A
  x.SA[0 : 1,1,1],       //设置初值为1,1,1 
  K# W$ B2 \! ^+ Z$ g
  i=netn[HFor("f"),x],   //拟牛顿法解方程 
. w, L+ \+ j0 _0 ^\\" E6 ?! \/ ]! d
  x.outm(),              //输出结果0 f+ D+ I8 E2 c8 q7 O. w. m

i                      //返回迭代次数/ S* {) H# Z8 i$ K* ^0 I8 P

}; 
2 }/ j9 d) {4 j, B' s2 {

结果：
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
& A* f4 @) K `8 x4 q4 w
clc
0 K* b9 z\" D3 }2 O7 Z5 U. U! @
tic
# \) U* U7 A. U7 n0 Q% G* n0 }
k = zeros(5,5); % //生成5×5全0矩阵
4 x& P U# s( x+ c( T2 n' |) K
% 循环计算以下程序段100000次：
; A( ~, N# k; N! M. a5 r+ i
for m = 1:100000
& m) R, Z* b& _' o3 r, _
a = rand(5,7);
\" n3 G! ~9 L$ Y$ @+ R
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化8 k q/ S2 V: _( T1 q
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);# G2 v* c- \* g L4 s3 o\" v
end+ |' y1 U- P0 V+ M: ~' O
k
: i+ S# ~8 t8 W' \0 U! w% [ ^$ a
toc
7 O7 Z d9 q }4 ?) }

复制代码

Forcal代码：

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

!using["math","sys"];
3 }* u! B( H+ @: E5 B: Y+ ~
mvar:1 u f0 M+ J! v' \$ v\\" R
t0=clock(),4 r( b9 U: r: [( ~
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0
+ J8 T5 B: I0 v' ~
i=0,(i<1000 00).while{ //循环计算1000 00次
& d* u/ v, V- c$ I% ?
oo{
& i$ G\\" u. }1 u\\" |+ R' I
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
7 r+ J$ S6 ~7 x7 P\\" U: U1 s
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)$ m( M, \. [5 S
},; F1 g9 z* d6 ^3 W4 F$ Q0 h6 @
i++/ v6 A7 U; J0 |
},, R, Z2 V. Y: Z3 p1 y2 \: {
k.outm(), //输出矩阵k，然后销毁k ~& G6 R \( r& e7 Z3 x4 r4 a
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"]; 
3 }* u! B( H+ @: E5 B: Y+ ~
mvar:1 u  f0 M+ J! v' \$ v\\" R

t0=clock(),4 r( b9 U: r: [( ~

oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0 
+ J8 T5 B: I0 v' ~
i=0,(i<1000 00).while{ //循环计算1000 00次 
& d* u/ v, V- c$ I% ?
  oo{ 
& i$ G\\" u. }1 u\\" |+ R' I
    a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化 
7 r+ J$ S6 ~7 x7 P\\" U: U1 s
    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)$ m( M, \. [5 S

},; F1 g9 z* d6 ^3 W4 F$ Q0 h6 @

i++/ v6 A7 U; J0 |

},, R, Z2 V. Y: Z3 p1 y2 \: {

k.outm(),             //输出矩阵k，然后销毁k  ~& G6 R  \( r& e7 Z3 x4 r4 a

[clock()-t0]/1000;    //得到计算时间，秒

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

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

!using["math","sys"];! _# q4 g5 ~1 w2 x, h$ P- B$ N\\" P4 K
(:t0,k,i,a,b)=
1 j0 Q% m7 K3 y( t1 U
{
$ B5 V' I# Q p8 i. t+ Z. @5 i/ B
t0=clock(),
. O! l# B+ s9 [( w: p' i+ Q
k=zeros[5,5],
- ^+ l7 Y3 F3 |6 c) D8 F# r
i=0,(i<1000 00).while{- O4 s/ k( t) V0 [7 h9 V; F8 b
oo{4 ?1 U1 D6 X2 ~4 b. A: c& ]\\" \
a=rand[5,7], b=rand[7,5],: _4 U2 H3 v% u4 ^\\" M
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
+ r6 `0 `5 ~0 W1 x# `
},1 j4 c: m8 T\\" u4 N
i++' D1 m. M3 K% t
},( h( l4 P/ c\\" r2 T* Z# d- I
k.outm().delete(),
: B9 G( l8 d# d- o( ]6 S
[clock()-t0]/1000) F# T' }0 x( v6 o; [
};

!using["math","sys"];! _# q4 g5 ~1 w2 x, h$ P- B$ N\\" P4 K

(:t0,k,i,a,b)= 
1 j0 Q% m7 K3 y( t1 U
{ 
$ B5 V' I# Q  p8 i. t+ Z. @5 i/ B
  t0=clock(), 
. O! l# B+ s9 [( w: p' i+ Q
  k=zeros[5,5], 
- ^+ l7 Y3 F3 |6 c) D8 F# r
  i=0,(i<1000 00).while{- O4 s/ k( t) V0 [7 h9 V; F8 b

oo{4 ?1 U1 D6 X2 ~4 b. A: c& ]\\" \

a=rand[5,7], b=rand[7,5],: _4 U2 H3 v% u4 ^\\" M

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
+ r6 `0 `5 ~0 W1 x# `
    },1 j4 c: m8 T\\" u4 N

i++' D1 m. M3 K% t

},( h( l4 P/ c\\" r2 T* Z# d- I

k.outm().delete(), 
: B9 G( l8 d# d- o( ]6 S
  [clock()-t0]/1000) F# T' }0 x( v6 o; [

};

在我的电脑上运行时间为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"];
5 N6 k3 @% @4 F; C
mvar:, U4 n% h; D7 q# R' B! I+ A
oo{
6 j0 J8 ^' e6 A5 n' A$ y
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],1 B6 N/ M2 q7 o1 ]5 [
a1.outm[5,1,1],. X( t2 O' Y, B! b. C* v
a2.outm[5,1,1],( B& H5 r, f: Z; ]8 G, l
a3.outm[5,1,1],* [* J; O' R9 b; F
a4.outm[5,1,1],
& D+ _+ P, {5 X4 N9 h2 _
a=a1+a2+a3+a4,% N6 `+ ?! p$ e T* u
a.outm[5,1,1],. g4 a! `: B6 o
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]. B3 n7 V5 s9 D$ v
};: J. n+ V' W% @ F, v/ D

!using["math"]; 
5 N6 k3 @% @4 F; C
mvar:, U4 n% h; D7 q# R' B! I+ A

oo{ 
6 j0 J8 ^' e6 A5 n' A$ y
  ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],1 B6 N/ M2 q7 o1 ]5 [

a1.outm[5,1,1],. X( t2 O' Y, B! b. C* v

a2.outm[5,1,1],( B& H5 r, f: Z; ]8 G, l

a3.outm[5,1,1],* [* J; O' R9 b; F

a4.outm[5,1,1], 
& D+ _+ P, {5 X4 N9 h2 _
  a=a1+a2+a3+a4,% N6 `+ ?! p$ e  T* u

a.outm[5,1,1],. g4 a! `: B6 o

Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]. B3 n7 V5 s9 D$ v

};: J. n+ V' W% @  F, v/ D

说明：
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"];
* ?4 }; g6 h: s0 O2 `* h
mvar:
/ P* B& v$ y2 Z$ Q
(:p1,p2,p3,a,b)=/ C7 Q- F3 C& M
{7 I7 `6 j1 ]4 N( W3 c3 x- e
oo{2 k8 x. w- g9 m) b# j: S; o7 w
a=array[1000].rand(),
+ j. `5 T9 k\\" e! B. F1 f
b=array[1000].rand(),
) j3 w6 l4 G& \/ _: Q
p1=array[1000,1000],' x2 L! B, C8 _& U2 J5 b\\" S! |* A
p2=array[1000,1000],
7 _; a4 k7 l/ q' i$ E5 e
p3=array[1000,1000],
/ u3 O E% P) c# T! _
t0=clock(),* O4 w+ W6 R5 a\\" ?
ndgrid(a,b,&A,&B),
- F6 a6 V& }9 F( I( t
p1.=A+B* t3 P+ u. S9 S3 H/ j N- u. o! L
},
3 X* L) ]) H+ ~- d
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},& u7 J) { I5 [& j& L/ _
lena=FCDLen(a),% y5 S; v- l; N1 a
lenb=FCDLen(b),4 }; f0 I! e2 m( F! n. ~- c
t0=clock(),. N% w& j i& l: z
m = lenb-1, (m>=0).while{
+ U2 Y7 v. m5 U$ D7 P
oo{p2(m,neg) = a+rn[b(m)]},
8 H# T' U1 {6 `* E
m--
8 _2 ?# U0 b# J5 f
},
$ i6 l9 `$ g7 o* J$ o0 u\\" Z
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},
# C$ p$ y6 R& f3 K. U; i1 F
t0=clock(),( K7 s8 \. `\\" E
m = lenb-1, (m>=0).while{
/ ?! g* j' p! g
n = lena-1, (n>=0).while{
6 C9 v, a7 _; e
//p3(m,n) = a(n)+b(m), //用这句还要慢一些# ]. I0 A! w4 K5 E* K
A(p3,m,n) = A(a,n)+A(b,m),- W* B+ J$ l* D& D* v
n--
: o6 B- I4 U\\" g; ^
},
. \9 v5 T( }\\" z1 \; x- i. W2 [) ~
m--
3 P, Q\\" C7 |5 D% v
},1 m: n5 Q* F C; s
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
3 F( p7 p1 ]6 d/ N
};* K0 y( l* q& f; ?0 Q

!using["math","sys"]; 
* ?4 }; g6 h: s0 O2 `* h
mvar: 
/ P* B& v$ y2 Z$ Q
(:p1,p2,p3,a,b)=/ C7 Q- F3 C& M

{7 I7 `6 j1 ]4 N( W3 c3 x- e

oo{2 k8 x. w- g9 m) b# j: S; o7 w

a=array[1000].rand(), 
+ j. `5 T9 k\\" e! B. F1 f
    b=array[1000].rand(), 
) j3 w6 l4 G& \/ _: Q
    p1=array[1000,1000],' x2 L! B, C8 _& U2 J5 b\\" S! |* A

p2=array[1000,1000], 
7 _; a4 k7 l/ q' i$ E5 e
    p3=array[1000,1000], 
/ u3 O  E% P) c# T! _
    t0=clock(),* O4 w+ W6 R5 a\\" ?

ndgrid(a,b,&A,&B), 
- F6 a6 V& }9 F( I( t
    p1.=A+B* t3 P+ u. S9 S3 H/ j  N- u. o! L

}, 
3 X* L) ]) H+ ~- d
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},& u7 J) {  I5 [& j& L/ _

lena=FCDLen(a),% y5 S; v- l; N1 a

lenb=FCDLen(b),4 }; f0 I! e2 m( F! n. ~- c

t0=clock(),. N% w& j  i& l: z

m = lenb-1, (m>=0).while{ 
+ U2 Y7 v. m5 U$ D7 P
    oo{p2(m,neg) = a+rn[b(m)]}, 
8 H# T' U1 {6 `* E
    m-- 
8 _2 ?# U0 b# J5 f
  }, 
$ i6 l9 `$ g7 o* J$ o0 u\\" Z
  printff{"\r\nfor1: {1,r}",[clock()-t0]/1000}, 
# C$ p$ y6 R& f3 K. U; i1 F
  t0=clock(),( K7 s8 \. `\\" E

m = lenb-1, (m>=0).while{ 
/ ?! g* j' p! g
    n = lena-1, (n>=0).while{ 
6 C9 v, a7 _; e
      //p3(m,n) = a(n)+b(m),  //用这句还要慢一些# ]. I0 A! w4 K5 E* K

A(p3,m,n) = A(a,n)+A(b,m),- W* B+ J$ l* D& D* v

n-- 
: o6 B- I4 U\\" g; ^
    }, 
. \9 v5 T( }\\" z1 \; x- i. W2 [) ~
    m-- 
3 P, Q\\" C7 |5 D% v
  },1 m: n5 Q* F  C; s

printff{"\r\nfor2: {1,r}",[clock()-t0]/1000} 
3 F( p7 p1 ]6 d/ N
};* K0 y( l* q& f; ?0 Q

结果：
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"];
]: C\\" n4 ?, S
mvar:
' ?7 Y3 b. e1 W8 D
t=clock(),! _# `; c& F8 K8 P: s
oo{
% X ^# w. }3 r
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
1 [: z% i1 \; o& o. j. l9 R
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]: j1 `8 b& @) y0 D
};; z* c. x2 w* C( d$ i0 k
[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)))));
) U+ G+ c8 X3 h\" f6 U1 k0 k4 j* z! f
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:
+ y. D' e! G5 W, p9 l
t=sys::clock();
@# k/ D* z9 b- ?+ [ r7 |
s=0,x=0, 9 Y. e3 T# J2 I& P' U6 @, y
while{x<=1, //while循环算法； ' ^' _& C, n6 E$ ]- Q
y=1, 4 S( B+ \& g/ d: B0 T8 d+ c% ?
while{y<=2,
, W. M3 j- R\" ~& b! G5 I
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), ( d; ?2 x# s# a% v% X
y=y+0.0009 , p$ ]0 `3 R z
}, \" x& }. r# f6 D
x=x+0.0011 ' y; b- X$ f% j& D4 E( G
},
3 L: j) _8 z( H) O1 c! c% i2 ^8 W
s;% c6 I. f- j( m5 a6 ]
[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 编辑

好深奥！~~~~
# U/ T+ }/ p0 X( [* yqbist 发表于 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"];* a, w+ n5 x\\" D+ M- G! _
(:a,b,k,t0)=
) |& _\\" Y3 W5 D3 p g
oo{% Y6 {% J, {# t( W; F
a=rand[1000,1000], b=rand[1000,1000],3 }2 A2 {/ a! C% l
t0=clock(),
2 R% ]& d% N& {. Q
k=a*b, //矩阵乘
' I' `& n7 ?4 G. q$ p0 M
k[1,3:5,9].outm()# f8 ^2 x1 k6 M9 l% F; B; B( _
},% U* I. w) a2 \9 K( z
[clock()-t0]/1000;
/ p: y$ ^# b! g7 [

结果：

238.447 247.837 247.065 248.105 247.058
$ d5 c$ r/ Z' X- r2 T
244.123 249.925 247.553 243.981 250.016
2 t4 y( I$ K7 O# }5 T
236.387 252.025 245.651 248.866 248.866 {\" I1 h) ^- ~: N
2.219 秒# _* P\" r( n; ?

复制代码

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

!using["math","sys","XSLSF"];
' C. O: D2 U# b0 G* w3 s
(:a,b,k,t0)=
; |# _3 @, d$ s4 ]+ @
oo{
3 G) I- ?$ E7 z% H
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],4 T/ C' p3 ]4 j6 n# q9 F- Y S( w
t0=clock(),) u# e4 m1 g+ p7 A
rmul[k:a,b], //矩阵乘
t+ k7 Z E* q E. K, L! N
k[1,3:5,9].outm()2 `/ N\\" p, ~2 I% x
},/ p' r4 N# f1 l1 L0 n8 }5 C
[clock()-t0]/1000;3 f7 x l+ K- ]6 V9 a. _

结果：

262.121 247.583 260.529 259.548 258.328, Q& O- X/ o `+ ^
255.413 246.563 254.356 250.548 251.509
4 Z, I3 W, U, | P' o4 H% M2 u1 G
256.152 247.725 259.444 250.827 249.816
9 H2 q% r# Q3 D3 P: Y
10.563 秒
! x% u0 I: p3 H( |

复制代码

我也说一句

发表

使用道具举报

返回列表

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

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

新人进步奖

QQ

电话咨询

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