查看: 5869|回复: 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"];3 c\\" U' u: O\\" h7 B* }. _
mvar:/ w+ Q! ^; w5 }7 F
oo{ //一般在oo函数中调用FcMath函数7 J4 H+ P' T\\" o% L0 x. z
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化; W: Z$ v/ E& P6 M& N\\" J
a.outm(), //输出矩阵a
# h8 |, I# p/ E/ o6 J- K, e6 S; S
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出4 A; O- B- J5 c: h: Y
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出8 w6 r# ^, q4 ~7 Y7 b; g
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出: b: C a* B- r2 Z0 b
};
# g+ b' `. y/ x# ~6 G

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.852615) p2 O\" B+ B7 o- J8 t# g0 G8 _ h
0.630646 0.927048 0.440308 0.162857 0.556854
' P& H& @\" ?# @ D+ L
0.43309 0.34552 0.563919 0.937164 0.2096415 q$ o i3 ~* J- I* ?5 {
0.603271 0.727676 0.130951 5.35736e-002 0.197937
* k! i4 _ {; y1 w\" m2 v5 |
0.576004 0.747589 1.17645e-002 0.363892 0.280777
( K8 K8 E( p4 M, N/ ?. {
0.646454 0.381088 0.58551 0.26387 0.936929 [( @8 e- ^% }8 G4 M& M
$ f& d( r& w3 ^. D$ e7 m
0.153259' p. D5 f4 X1 i/ h/ t/ E: Z\" R
0.1628572 Z\" b4 _$ K, a; v, K
0.937164
6 Z$ X5 o7 H6 C1 ?
5.35736e-002% j; e, { {4 R' C* s) [* ~
0.363892; ?6 ~3 t\" s- R
0.263879 T7 ^0 `2 ?6 B5 l: X/ ~
) o4 f) {* E+ z, {
0.603271 0.727676 0.130951 5.35736e-002 0.197937
/ h, H3 U8 E0 P0 [
\" Q7 Z& U1 z: V6 h& D; i! g
0.130951 5.35736e-002; |! w7 y9 f I; {3 f2 ~# W
1.17645e-002 0.363892
P7 J! G* }; z9 o9 H
0.58551 0.26387
: w* X! ~, f# ~; T
( ]) e6 h: m# L2 c

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义
8 f6 S8 R- v5 H' S9 I+ Z p. q
{
1 y\\" N1 q7 @+ k- d; G! E- l( j
y1=x1*x1+x2*x2+x3*x3-1.0,& v5 ^% }7 Y2 c a, v# E' E
y2=2.0*x1*x1+x2*x2-4.0*x3,
\\" s. L# r% [* Q* Q
y3=3.0*x1*x1-4.0*x2+x3*x3
. X; T. x\\" n4 ~ i+ z8 P
};
* n% X, z& r+ W7 n7 ]. t5 \
!using["math","sys"];
3 o& x/ O' j, u- G
mvar:$ m$ _' \2 @: W. F7 M. n G Y
oo{
: `5 J T( B; |, Y1 L3 n
x=array(3),! K\\" c3 d9 y1 |4 c3 ^
x.SA[0 : 1,1,1], //设置初值为1,1,1
% I: e, N- `% k# m [
i=netn[HFor("f"),x], //拟牛顿法解方程
6 n$ e- D8 O3 |, U
x.outm(), //输出结果
* ^3 S4 L\\" h2 h
i //返回迭代次数2 o# S% T! o9 r$ ~! _2 N/ e3 U
};
: P* B7 g# A0 h! g$ g

f(x1,x2,x3,y1,y2,y3)=      //函数定义 
8 f6 S8 R- v5 H' S9 I+ Z  p. q
{ 
1 y\\" N1 q7 @+ k- d; G! E- l( j
    y1=x1*x1+x2*x2+x3*x3-1.0,& v5 ^% }7 Y2 c  a, v# E' E

y2=2.0*x1*x1+x2*x2-4.0*x3, 
\\" s. L# r% [* Q* Q
    y3=3.0*x1*x1-4.0*x2+x3*x3 
. X; T. x\\" n4 ~  i+ z8 P
}; 
* n% X, z& r+ W7 n7 ]. t5 \
!using["math","sys"]; 
3 o& x/ O' j, u- G
mvar:$ m$ _' \2 @: W. F7 M. n  G  Y

oo{ 
: `5 J  T( B; |, Y1 L3 n
  x=array(3),! K\\" c3 d9 y1 |4 c3 ^

x.SA[0 : 1,1,1],       //设置初值为1,1,1 
% I: e, N- `% k# m  [
  i=netn[HFor("f"),x],   //拟牛顿法解方程 
6 n$ e- D8 O3 |, U
  x.outm(),              //输出结果 
* ^3 S4 L\\" h2 h
  i                      //返回迭代次数2 o# S% T! o9 r$ ~! _2 N/ e3 U

}; 
: P* B7 g# A0 h! g$ g

结果：
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
, i/ P* O$ `6 A ~3 b1 h
clc
) v' u; R# V$ l# N\" H! `
tic& k\" w. D4 z9 h
k = zeros(5,5); % //生成5×5全0矩阵' z! O3 i7 p! P7 M% p; z2 ?
% 循环计算以下程序段100000次：
! ^+ ?: q- i6 c8 y' Z# _
for m = 1:100000
2 o% }! S9 d\" o! \$ A5 R2 `
a = rand(5,7);
8 t ?* V: G\" M! U9 G/ U1 ?2 f
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
* j) h. o! A$ w/ n$ P
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);- a! T- J3 d5 g B: w
end4 w\" P1 W4 E0 g6 ^ r
k4 N/ B* u2 }: R/ J0 [9 f
toc
7 n! k+ Q7 }2 V) c9 \

复制代码

Forcal代码：

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

!using["math","sys"];
5 |' X* [5 J9 H1 R s) ~8 L5 `
mvar:
$ K( b- P' t: ^6 t7 F! J\\" T6 R
t0=clock(),9 ^ c* O) O$ I( ]; z. f
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0
$ i6 g, @, q. z9 f8 A5 @8 W
i=0,(i<1000 00).while{ //循环计算1000 00次
9 `; ? h7 {: P( S
oo{
8 j+ D; n0 N1 U6 C
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
. q8 }- z/ H( H\\" R
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)
& J4 Q2 P# S6 E* `* m; f) J
},& x, w0 Q2 S, _: W
i++- m3 P8 X3 V' m+ Y9 A0 e- E9 W; j& h
},2 w8 ^3 _: }; u+ E' B- _
k.outm(), //输出矩阵k，然后销毁k3 F, S7 g/ t: @1 c7 m1 T
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"]; 
5 |' X* [5 J9 H1 R  s) ~8 L5 `
mvar: 
$ K( b- P' t: ^6 t7 F! J\\" T6 R
t0=clock(),9 ^  c* O) O$ I( ]; z. f

oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0 
$ i6 g, @, q. z9 f8 A5 @8 W
i=0,(i<1000 00).while{ //循环计算1000 00次 
9 `; ?  h7 {: P( S
  oo{ 
8 j+ D; n0 N1 U6 C
    a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化 
. q8 }- z/ H( H\\" R
    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) 
& J4 Q2 P# S6 E* `* m; f) J
  },& x, w0 Q2 S, _: W

i++- m3 P8 X3 V' m+ Y9 A0 e- E9 W; j& h

},2 w8 ^3 _: }; u+ E' B- _

k.outm(),             //输出矩阵k，然后销毁k3 F, S7 g/ t: @1 c7 m1 T

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

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

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

!using["math","sys"];7 I6 I7 E8 t* }3 {: y
(:t0,k,i,a,b)=
: Q0 N! }% c% V. ]. ?) X% S
{
0 m7 @8 ~\\" u+ K5 T
t0=clock(),
7 S# R0 f. c$ a0 z4 l
k=zeros[5,5],\\" b3 n6 l1 U! l\\" }
i=0,(i<1000 00).while{' w5 d0 c6 g: ?/ U5 G- }
oo{8 R, r9 f5 D3 V3 n2 e5 M( d
a=rand[5,7], b=rand[7,5],1 J) F5 |6 P\\" w: }, T! m4 y
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
\\" q\\" U0 x3 H# V\\" W: z) x
},, U4 o/ q# S\\" J; w% z
i++
/ |7 e7 K. z. _! ~' J& T: k
},
+ p/ u0 j @, D. L
k.outm().delete(),
0 y% {2 ]% I# h0 B
[clock()-t0]/1000$ \; ]9 g& t% _& s
};

!using["math","sys"];7 I6 I7 E8 t* }3 {: y

(:t0,k,i,a,b)= 
: Q0 N! }% c% V. ]. ?) X% S
{ 
0 m7 @8 ~\\" u+ K5 T
  t0=clock(), 
7 S# R0 f. c$ a0 z4 l
  k=zeros[5,5],\\" b3 n6 l1 U! l\\" }

i=0,(i<1000 00).while{' w5 d0 c6 g: ?/ U5 G- }

oo{8 R, r9 f5 D3 V3 n2 e5 M( d

a=rand[5,7], b=rand[7,5],1 J) F5 |6 P\\" w: }, T! m4 y

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
\\" q\\" U0 x3 H# V\\" W: z) x
    },, U4 o/ q# S\\" J; w% z

i++ 
/ |7 e7 K. z. _! ~' J& T: k
  }, 
+ p/ u0 j  @, D. L
  k.outm().delete(), 
0 y% {2 ]% I# h0 B
  [clock()-t0]/1000$ \; ]9 g& t% _& s

};

在我的电脑上运行时间为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"];
# s9 U: b! i* F0 y4 A) ~5 Z
mvar:
7 E9 G' ?7 j9 o: E! F9 x\\" T
oo{
Z7 G% Q5 R0 i9 ?! c3 M
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
. ^6 ?1 u# @' q
a1.outm[5,1,1],
3 ] b$ @* L* | D
a2.outm[5,1,1],
& ~+ I$ k# Y! `) b+ g' |
a3.outm[5,1,1],
$ W: y# ~3 y6 p$ M$ j
a4.outm[5,1,1],1 B& W1 x$ ~& X/ R( ?
a=a1+a2+a3+a4,
0 ?3 K; A. c7 @- M% E; q
a.outm[5,1,1],
! P' c0 C' t, ~1 l. ]4 T1 p
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]
* K; f) t( D5 Q0 c/ }6 f |9 c: K% k
};% o t* D' w! x+ b5 s

!using["math"]; 
# s9 U: b! i* F0 y4 A) ~5 Z
mvar: 
7 E9 G' ?7 j9 o: E! F9 x\\" T
oo{ 
  Z7 G% Q5 R0 i9 ?! c3 M
  ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4], 
. ^6 ?1 u# @' q
  a1.outm[5,1,1], 
3 ]  b$ @* L* |  D
  a2.outm[5,1,1], 
& ~+ I$ k# Y! `) b+ g' |
  a3.outm[5,1,1], 
$ W: y# ~3 y6 p$ M$ j
  a4.outm[5,1,1],1 B& W1 x$ ~& X/ R( ?

a=a1+a2+a3+a4, 
0 ?3 K; A. c7 @- M% E; q
  a.outm[5,1,1], 
! P' c0 C' t, ~1 l. ]4 T1 p
  Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[] 
* K; f) t( D5 Q0 c/ }6 f  |9 c: K% k
};% o  t* D' w! x+ b5 s

说明：
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"];8 [\\" b7 S$ N8 U& F2 P
mvar:
) ]6 M4 Z2 Y# J( r W% F. ]
(:p1,p2,p3,a,b)=
; p3 u; @: d$ U' }\\" I
{
1 `/ D g! C( t
oo{
7 d0 b4 Q' k6 ?# A( }
a=array[1000].rand(),
5 ?2 N# g; a4 V: l8 D. q& |
b=array[1000].rand(),
3 w+ G) y. O4 I% T: C' x/ S
p1=array[1000,1000],) u5 n4 X3 Z* s1 V\\" g' f- @/ U/ v
p2=array[1000,1000],
2 _2 a* @( V1 h/ N: X
p3=array[1000,1000],2 q5 `$ p\\" `: D4 n
t0=clock(),. _/ R6 v; W* x, m- K; u
ndgrid(a,b,&A,&B),
- \; K4 g! e( b9 h; \% K7 Y
p1.=A+B6 n1 s7 N m9 Q1 r6 q- J o! W
},
- K0 T* y& s/ U& I/ [
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},
d0 X9 k. P: q) r/ I7 E\\" P4 I
lena=FCDLen(a),, D1 B8 C; Q9 n* o6 _5 d6 e% m. x
lenb=FCDLen(b),
/ H# l6 [) N6 Y. @
t0=clock(),
& K% f m6 u+ N% Z
m = lenb-1, (m>=0).while{
% N8 ]1 j( l& h& V% C
oo{p2(m,neg) = a+rn[b(m)]},! _\\" Q8 k% q& I7 \ ^
m--
2 [ Y8 C3 p/ y7 y
},
- J. |6 S+ u, h
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},
! N1 ~8 A, `( Z r/ z
t0=clock(),2 t\\" s- p' `' f% [0 h) O' y, k
m = lenb-1, (m>=0).while{
8 s2 N$ b5 g- h: A/ c) f
n = lena-1, (n>=0).while{8 @7 d5 t& P7 g: k
//p3(m,n) = a(n)+b(m), //用这句还要慢一些* _& T0 {( C4 S
A(p3,m,n) = A(a,n)+A(b,m),
. E- \* X\\" l, B3 S
n--; b5 O5 N Q7 Q1 R
},
7 b( N# p8 ^* i! J2 E, a5 F A
m--8 S( f2 g, S0 X7 f$ u
},- y1 I' A* q; e& q# H. E' b( b
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
};% W% m9 t- X, [/ F/ X: E6 \- M0 Z+ Y; _: M

!using["math","sys"];8 [\\" b7 S$ N8 U& F2 P

mvar: 
) ]6 M4 Z2 Y# J( r  W% F. ]
(:p1,p2,p3,a,b)= 
; p3 u; @: d$ U' }\\" I
{ 
1 `/ D  g! C( t
  oo{ 
7 d0 b4 Q' k6 ?# A( }
    a=array[1000].rand(), 
5 ?2 N# g; a4 V: l8 D. q& |
    b=array[1000].rand(), 
3 w+ G) y. O4 I% T: C' x/ S
    p1=array[1000,1000],) u5 n4 X3 Z* s1 V\\" g' f- @/ U/ v

p2=array[1000,1000], 
2 _2 a* @( V1 h/ N: X
    p3=array[1000,1000],2 q5 `$ p\\" `: D4 n

t0=clock(),. _/ R6 v; W* x, m- K; u

ndgrid(a,b,&A,&B), 
- \; K4 g! e( b9 h; \% K7 Y
    p1.=A+B6 n1 s7 N  m9 Q1 r6 q- J  o! W

}, 
- K0 T* y& s/ U& I/ [
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000}, 
  d0 X9 k. P: q) r/ I7 E\\" P4 I
  lena=FCDLen(a),, D1 B8 C; Q9 n* o6 _5 d6 e% m. x

lenb=FCDLen(b), 
/ H# l6 [) N6 Y. @
  t0=clock(), 
& K% f  m6 u+ N% Z
  m = lenb-1, (m>=0).while{ 
% N8 ]1 j( l& h& V% C
    oo{p2(m,neg) = a+rn[b(m)]},! _\\" Q8 k% q& I7 \  ^

m-- 
2 [  Y8 C3 p/ y7 y
  }, 
- J. |6 S+ u, h
  printff{"\r\nfor1: {1,r}",[clock()-t0]/1000}, 
! N1 ~8 A, `( Z  r/ z
  t0=clock(),2 t\\" s- p' `' f% [0 h) O' y, k

m = lenb-1, (m>=0).while{ 
8 s2 N$ b5 g- h: A/ c) f
    n = lena-1, (n>=0).while{8 @7 d5 t& P7 g: k

//p3(m,n) = a(n)+b(m),  //用这句还要慢一些* _& T0 {( C4 S

A(p3,m,n) = A(a,n)+A(b,m), 
. E- \* X\\" l, B3 S
      n--; b5 O5 N  Q7 Q1 R

}, 
7 b( N# p8 ^* i! J2 E, a5 F  A
    m--8 S( f2 g, S0 X7 f$ u

},- y1 I' A* q; e& q# H. E' b( b

printff{"\r\nfor2: {1,r}",[clock()-t0]/1000} 
! q9 @( E5 s% g) u\\" b2 E* \1 s6 H7 R% f
};% W% m9 t- X, [/ F/ X: E6 \- M0 Z+ Y; _: M

结果：
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"];
( E0 y) L; [% y8 j- b9 t1 k+ I
mvar:. Z0 U1 S7 X8 X Y
t=clock(),# s1 P- F( m/ k k
oo{( Q5 S# D0 k6 u/ p8 o
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],6 ?& S4 i' p# E( G) q: |
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]% [5 W! e5 j, T2 O' | ]& q+ y+ l; s
};7 O7 M( }5 O+ z7 O$ T
[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))))); ) _# ]# R$ c5 G/ [' d: J) I
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:
4 i& w& u0 g( ?! S [( ?$ J
t=sys::clock();$ C\" {6 C. ~( _1 a) Q4 C
s=0,x=0, ; M, e7 o1 m' i3 l% M
while{x<=1, //while循环算法；
\" Y5 h$ t/ P8 L# J' c0 E. t; K A% S
y=1, ' o! o8 D- ?% ?- k+ X
while{y<=2,
6 v1 y$ j' Z c$ d7 c h; s
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
4 ~\" H) P b D/ b4 S! `5 B( ]
y=y+0.0009
5 z( y\" ^6 p. c4 ?5 I+ f, V L: o& m
}, ) S8 q3 B- @+ ?1 R1 x
x=x+0.0011 3 A2 o6 D% z Z9 C- d6 k. X
},
$ X: T/ o$ O$ W9 c8 L
s;+ V5 }( ?' z; N9 ?& l& M! G, I
[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 编辑

好深奥！~~~~2 Y5 R# C u8 U. g
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"];/ w. R5 w+ @2 s* D/ [4 z. ]
(:a,b,k,t0)=9 a* L+ _7 E\\" H& ?
oo{( q& }( |1 b2 F7 T
a=rand[1000,1000], b=rand[1000,1000],
* \8 q S2 s$ c0 i
t0=clock(),
\\" N- U% K$ B3 l6 Y6 M9 m/ L' R
k=a*b, //矩阵乘
5 G) d' Q7 t9 N) ]
k[1,3:5,9].outm()7 [7 [* W2 E$ H) c
},0 K, N1 }6 x; T6 H1 w
[clock()-t0]/1000;6 e/ ?1 m% _6 s. N8 J\\" q

结果：

238.447 247.837 247.065 248.105 247.058
244.123 249.925 247.553 243.981 250.016' j1 F8 H% X3 P$ r+ _: n3 [, j
236.387 252.025 245.651 248.866 248.866
9 i8 Y9 A6 o+ F4 z, M1 X
2.219 秒9 T9 h1 e$ G; ?0 M+ g

复制代码

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

!using["math","sys","XSLSF"];' a, @! U' b, I1 b3 j: w( v( ?, D
(:a,b,k,t0)=- ~3 c9 d$ F& r: ]% j
oo{# s& l$ u/ n+ ~! N
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],/ R9 J9 V0 e2 I1 W- ^2 O! ]
t0=clock(),% l. j7 P3 V: L% k+ B
rmul[k:a,b], //矩阵乘
) q/ s/ @\\" c8 z3 r2 `& ?9 A
k[1,3:5,9].outm()* y9 `0 D0 Y; {4 ^
},. m7 T! G0 Y' @$ J; Q% p( t\\" e* p
[clock()-t0]/1000;
3 a, p8 o; y* g% q0 a0 w! K

结果：

262.121 247.583 260.529 259.548 258.328
9 s2 ?\" M) S; M7 N0 N. T! ~8 a
255.413 246.563 254.356 250.548 251.509+ N6 K# [4 i% n+ Z* ]8 Q) _
256.152 247.725 259.444 250.827 249.816
. G% f+ ?0 x$ ]+ q d! m. R
10.563 秒
6 I7 L: l r9 i6 C\" w' ] y

复制代码

我也说一句

发表

使用道具举报

返回列表

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

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

新人进步奖

QQ

电话咨询

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