查看: 5604|回复: 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"];9 E0 o9 C! M\\" M# m
mvar:, a! h- Y* r* d! h1 m# L
oo{ //一般在oo函数中调用FcMath函数* k' D- t: j\\" j4 E) [; _- O\\" n( I
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化
, z! t, W: I* f9 W9 `2 M
a.outm(), //输出矩阵a
, G0 Z( }0 r* X
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出
) ?; O' X4 ?& W
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出
U% C4 A9 C- |
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出' U! u\\" z6 s0 M9 S0 m* J( Z+ G7 Z8 d P
};5 A. I\\" }! D: M! o; q

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.8526150 d- ?) ?/ m0 a6 u( y. ^/ e% ?
0.630646 0.927048 0.440308 0.162857 0.5568540 Q8 w- r( L1 L9 P( T6 R9 O' e2 f
0.43309 0.34552 0.563919 0.937164 0.2096413 ^, Z- T/ X1 D: r
0.603271 0.727676 0.130951 5.35736e-002 0.1979371 m' i\" U0 S# {. A0 g
0.576004 0.747589 1.17645e-002 0.363892 0.280777
/ O8 J9 S6 J3 k! u\" U\" I
0.646454 0.381088 0.58551 0.26387 0.93692
: y2 f! Q# R% u3 H, q) ?. Z
\" |- P- K. U* A& v
0.153259
- J) x; V' M, t3 m' C7 P. Z2 _
0.162857
- s4 H/ |. w' R
0.937164' a( |* g/ y4 _% R5 X$ m* Y
5.35736e-002
2 O7 s/ Z1 E: F0 ^1 Q
0.363892
. i0 R8 |$ s5 i5 Z8 ^
0.26387
* Y1 o% L; S) c. C% X8 Z
0 ^/ B1 m+ ~$ u( e( u3 T, A( k5 s
0.603271 0.727676 0.130951 5.35736e-002 0.197937* W$ q9 X' U3 t) W) m7 r# C
2 j* p+ V# P, w7 K
0.130951 5.35736e-002/ A3 z* l- b A, z7 e. Q% e6 e
1.17645e-002 0.363892
- P* Y$ D: `# ^' |; j\" }
0.58551 0.26387
* v- }' i$ b( ]1 H
4 S1 `9 F8 J1 T- g* ~+ @

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义
\\" n% c' o! K8 }\\" W7 e
{
/ l% g V7 q9 h; _
y1=x1*x1+x2*x2+x3*x3-1.0,
7 E( c/ r8 R\\" v+ H7 D+ \) D) @
y2=2.0*x1*x1+x2*x2-4.0*x3,5 R+ z; w\\" R. y( v3 V k
y3=3.0*x1*x1-4.0*x2+x3*x3* Q2 p( T% v& ~& |1 _' C1 R8 v& d
};
7 E W. h2 K1 H- V8 C$ S
!using["math","sys"];9 K% J* q4 k/ R0 C
mvar:0 x. C3 p+ |7 o2 g: Z' w
oo{
g* ~/ P1 k1 U: q1 e5 i2 _
x=array(3),
u& l# `7 U2 V/ L
x.SA[0 : 1,1,1], //设置初值为1,1,1$ Z6 m; h1 T$ E8 [8 f
i=netn[HFor("f"),x], //拟牛顿法解方程- l7 G, G N% ?& z) {5 }! Y7 ^
x.outm(), //输出结果; q' Z& Y8 I L. K: k! m
i //返回迭代次数
6 H* D2 N\\" k9 x- ~0 h7 b
};) `' j2 y' b. @\\" f/ ?\\" t

f(x1,x2,x3,y1,y2,y3)=      //函数定义 
\\" n% c' o! K8 }\\" W7 e
{ 
/ l% g  V7 q9 h; _
    y1=x1*x1+x2*x2+x3*x3-1.0, 
7 E( c/ r8 R\\" v+ H7 D+ \) D) @
    y2=2.0*x1*x1+x2*x2-4.0*x3,5 R+ z; w\\" R. y( v3 V  k

y3=3.0*x1*x1-4.0*x2+x3*x3* Q2 p( T% v& ~& |1 _' C1 R8 v& d

}; 
7 E  W. h2 K1 H- V8 C$ S
!using["math","sys"];9 K% J* q4 k/ R0 C

mvar:0 x. C3 p+ |7 o2 g: Z' w

oo{ 
  g* ~/ P1 k1 U: q1 e5 i2 _
  x=array(3), 
  u& l# `7 U2 V/ L
  x.SA[0 : 1,1,1],       //设置初值为1,1,1$ Z6 m; h1 T$ E8 [8 f

i=netn[HFor("f"),x],   //拟牛顿法解方程- l7 G, G  N% ?& z) {5 }! Y7 ^

x.outm(),              //输出结果; q' Z& Y8 I  L. K: k! m

i                      //返回迭代次数 
6 H* D2 N\\" k9 x- ~0 h7 b
};) `' j2 y' b. @\\" f/ ?\\" t

结果：
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* Q! J, `# I. X) `
clc0 }# E: u# @3 H6 W( Z
tic
4 _0 }, y\" \$ Z0 f% E7 }; {1 K9 |
k = zeros(5,5); % //生成5×5全0矩阵
, E, Z\" b* ^0 H8 @3 O9 [
% 循环计算以下程序段100000次：
% j( T; l) e' s& U& L% w
for m = 1:1000002 }2 q9 F; I3 \$ L X9 w
a = rand(5,7);# ]8 @* k% _2 Y( s
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
\" v* b3 V% K$ c7 m
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);9 o( o1 s+ s7 J5 k\" ^; z
end
: J3 r v2 c* t* J; C- q( _
k% ~\" M* R& Y5 _\" Y A
toc
8 b0 U8 ]3 o9 l, e+ y

复制代码

Forcal代码：

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

!using["math","sys"];1 o3 J+ [9 F5 c c8 I
mvar:
! t+ A3 o1 X. j/ ^$ b1 g
t0=clock(),1 W. u' o/ i/ s7 L- @# D0 A2 w
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0
6 \- J4 D7 H7 c4 q& j\\" d
i=0,(i<1000 00).while{ //循环计算1000 00次
, Q. w4 K- r+ V
oo{
8 b8 C0 j) s: Y8 q$ W
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
. T# q* @: ^0 [\\" a0 Q0 Z$ T1 j
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)
9 I( ~0 i* ^% ~' L, @
},
; n: H x% W& d: I; X; N0 b
i++/ c8 h3 Q3 b) o* i+ }5 ?( `: Q5 e. a
},
6 B8 k1 Q* V* I% N9 q2 i
k.outm(), //输出矩阵k，然后销毁k
3 r: `* \: u# |) ~, u. R
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"];1 o3 J+ [9 F5 c  c8 I

mvar: 
! t+ A3 o1 X. j/ ^$ b1 g
t0=clock(),1 W. u' o/ i/ s7 L- @# D0 A2 w

oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0 
6 \- J4 D7 H7 c4 q& j\\" d
i=0,(i<1000 00).while{ //循环计算1000 00次 
, Q. w4 K- r+ V
  oo{ 
8 b8 C0 j) s: Y8 q$ W
    a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化 
. T# q* @: ^0 [\\" a0 Q0 Z$ T1 j
    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) 
9 I( ~0 i* ^% ~' L, @
  }, 
; n: H  x% W& d: I; X; N0 b
  i++/ c8 h3 Q3 b) o* i+ }5 ?( `: Q5 e. a

}, 
6 B8 k1 Q* V* I% N9 q2 i
k.outm(),             //输出矩阵k，然后销毁k 
3 r: `* \: u# |) ~, u. R
[clock()-t0]/1000;    //得到计算时间，秒

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

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

!using["math","sys"];, V5 e$ c\\" \\\" K. m! o
(:t0,k,i,a,b)=6 g ^5 \& W. k, B6 @/ L
{
/ q0 @% L( F2 h8 ]1 F8 f. q
t0=clock(),
9 D5 l# \% ~$ d7 [! E
k=zeros[5,5],
: L3 b) o, i6 K9 T( t% X
i=0,(i<1000 00).while{7 E* N7 ~+ }! K9 Y: A\\" Z
oo{% J0 G2 { c4 Z
a=rand[5,7], b=rand[7,5],
/ }6 K/ |6 z( A8 I. F+ I( C5 m
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) Q# V- T \9 F+ }7 }! B$ C( g\\" l
},1 R/ F: h3 x5 x. Y
i++
. t! V/ t. {+ }' g; L0 y( s
},; G, j( S# J0 H7 ?! i
k.outm().delete(),: R: W\\" [% r6 Z# P6 v
[clock()-t0]/1000, E7 U: W4 c2 `( @8 W\\" e
};

!using["math","sys"];, V5 e$ c\\" \\\" K. m! o

(:t0,k,i,a,b)=6 g  ^5 \& W. k, B6 @/ L

{ 
/ q0 @% L( F2 h8 ]1 F8 f. q
  t0=clock(), 
9 D5 l# \% ~$ d7 [! E
  k=zeros[5,5], 
: L3 b) o, i6 K9 T( t% X
  i=0,(i<1000 00).while{7 E* N7 ~+ }! K9 Y: A\\" Z

oo{% J0 G2 {  c4 Z

a=rand[5,7], b=rand[7,5], 
/ }6 K/ |6 z( A8 I. F+ I( C5 m
      k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)  Q# V- T  \9 F+ }7 }! B$ C( g\\" l

},1 R/ F: h3 x5 x. Y

i++ 
. t! V/ t. {+ }' g; L0 y( s
  },; G, j( S# J0 H7 ?! i

k.outm().delete(),: R: W\\" [% r6 Z# P6 v

[clock()-t0]/1000, E7 U: W4 c2 `( @8 W\\" e

};

在我的电脑上运行时间为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"];' c: v, H. x. k) C8 R8 d
mvar:
2 G, P/ s# _5 g9 s9 C# S# x
oo{- s; p, S+ l- \; \9 u: z
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],
7 w y8 e5 q% M4 K) H* U* ^
a1.outm[5,1,1],
# p. [# e/ C0 P# }- p& w& e8 p
a2.outm[5,1,1],5 F; I& N' V- R# j! ?/ Z
a3.outm[5,1,1],1 o! y3 @) b: \8 q/ D, A8 n
a4.outm[5,1,1],. g' K$ j. I1 u
a=a1+a2+a3+a4,5 f' R2 S+ }\\" x# h\\" O\\" E9 f4 ^4 v
a.outm[5,1,1],
. ?( Z% G4 f! I! k7 a8 I Y- i
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]% r7 r, H0 G4 q- K
};

!using["math"];' c: v, H. x. k) C8 R8 d

mvar: 
2 G, P/ s# _5 g9 s9 C# S# x
oo{- s; p, S+ l- \; \9 u: z

ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4], 
7 w  y8 e5 q% M4 K) H* U* ^
  a1.outm[5,1,1], 
# p. [# e/ C0 P# }- p& w& e8 p
  a2.outm[5,1,1],5 F; I& N' V- R# j! ?/ Z

a3.outm[5,1,1],1 o! y3 @) b: \8 q/ D, A8 n

a4.outm[5,1,1],. g' K$ j. I1 u

a=a1+a2+a3+a4,5 f' R2 S+ }\\" x# h\\" O\\" E9 f4 ^4 v

a.outm[5,1,1], 
. ?( Z% G4 f! I! k7 a8 I  Y- i
  Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]% r7 r, H0 G4 q- K

}; 
3 b% n% X1 y$ N  Z\\" I9 h. Q2 _/ R! X5 N

说明：
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 m' @* t# F! _$ d
mvar:
9 m3 q6 }) U, V0 O4 k9 n! w
(:p1,p2,p3,a,b)=' C: v6 u7 j7 F. w- V
{
, @, e- s3 n4 a4 n. [
oo{3 e$ ~9 A) |: d\\" A E Q; q
a=array[1000].rand(),0 |# d3 f! `\\" H
b=array[1000].rand(),3 w; I# B9 X- _6 w# \
p1=array[1000,1000],( W5 c) a7 ?$ ~! T( r/ v* J
p2=array[1000,1000],3 F3 [( ?5 K# o% u( V
p3=array[1000,1000],; A5 `0 f' T7 K$ [
t0=clock(),. v5 J8 S( J( s5 ]& i% G
ndgrid(a,b,&A,&B),( P0 t9 L; c' x: @
p1.=A+B i c' v ^' Y2 w/ V A
},
. l, ]3 v/ J( a) O) t1 \1 I
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},
, o\\" N! e5 A' w0 V\\" j
lena=FCDLen(a),
, O* e4 }$ F0 W) ]8 m: u# u% C
lenb=FCDLen(b),; ^& h4 A g# m) i6 B3 r& G
t0=clock(),
* y: E+ Q4 a( f# \' `+ ]
m = lenb-1, (m>=0).while{! b7 r S B5 J
oo{p2(m,neg) = a+rn[b(m)]},
# A. C! n8 ^6 f$ p\\" h8 t
m--
9 j! J0 Y# @7 B. b7 `5 a0 [; M
},
3 o5 l$ {( d2 E- V
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},/ {: K+ V7 j$ q+ j& ~4 u! X
t0=clock(),2 r6 r. `1 J\\" s! f) h1 C
m = lenb-1, (m>=0).while{4 d2 _# O0 c' G9 J
n = lena-1, (n>=0).while{
5 }, ?- ?\\" o0 }) ?9 v3 Y. g. Q& v4 v
//p3(m,n) = a(n)+b(m), //用这句还要慢一些
2 W3 q' O, v5 K5 d# e) Y0 o
A(p3,m,n) = A(a,n)+A(b,m),* j; f6 L1 o- o: c
n--
X- ^\\" l% x6 W- V
},
d& k: n* v3 u% {+ w1 [- r
m--) U# t o; N k! Q
},
; A8 ?8 K: q: N% E5 i6 s
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
* m' [) ~1 X\\" _0 i, ?! g
};1 s! o( D9 y1 h( Z. y2 L9 x

!using["math","sys"];5 m' @* t# F! _$ d

mvar: 
9 m3 q6 }) U, V0 O4 k9 n! w
(:p1,p2,p3,a,b)=' C: v6 u7 j7 F. w- V

{ 
, @, e- s3 n4 a4 n. [
  oo{3 e$ ~9 A) |: d\\" A  E  Q; q

a=array[1000].rand(),0 |# d3 f! `\\" H

b=array[1000].rand(),3 w; I# B9 X- _6 w# \

p1=array[1000,1000],( W5 c) a7 ?$ ~! T( r/ v* J

p2=array[1000,1000],3 F3 [( ?5 K# o% u( V

p3=array[1000,1000],; A5 `0 f' T7 K$ [

t0=clock(),. v5 J8 S( J( s5 ]& i% G

ndgrid(a,b,&A,&B),( P0 t9 L; c' x: @

p1.=A+B  i  c' v  ^' Y2 w/ V  A

}, 
. l, ]3 v/ J( a) O) t1 \1 I
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000}, 
, o\\" N! e5 A' w0 V\\" j
  lena=FCDLen(a), 
, O* e4 }$ F0 W) ]8 m: u# u% C
  lenb=FCDLen(b),; ^& h4 A  g# m) i6 B3 r& G

t0=clock(), 
* y: E+ Q4 a( f# \' `+ ]
  m = lenb-1, (m>=0).while{! b7 r  S  B5 J

oo{p2(m,neg) = a+rn[b(m)]}, 
# A. C! n8 ^6 f$ p\\" h8 t
    m-- 
9 j! J0 Y# @7 B. b7 `5 a0 [; M
  }, 
3 o5 l$ {( d2 E- V
  printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},/ {: K+ V7 j$ q+ j& ~4 u! X

t0=clock(),2 r6 r. `1 J\\" s! f) h1 C

m = lenb-1, (m>=0).while{4 d2 _# O0 c' G9 J

n = lena-1, (n>=0).while{ 
5 }, ?- ?\\" o0 }) ?9 v3 Y. g. Q& v4 v
      //p3(m,n) = a(n)+b(m),  //用这句还要慢一些 
2 W3 q' O, v5 K5 d# e) Y0 o
      A(p3,m,n) = A(a,n)+A(b,m),* j; f6 L1 o- o: c

n-- 
  X- ^\\" l% x6 W- V
    }, 
  d& k: n* v3 u% {+ w1 [- r
    m--) U# t  o; N  k! Q

}, 
; A8 ?8 K: q: N% E5 i6 s
  printff{"\r\nfor2: {1,r}",[clock()-t0]/1000} 
* m' [) ~1 X\\" _0 i, ?! g
};1 s! o( D9 y1 h( Z. y2 L9 x

结果：
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"];1 Y, I, W _4 f) j( y8 u
mvar:8 o; k% J3 D5 {: z# ^, T: C% t
t=clock(),$ v$ R' U' d\\" z) N0 U& x
oo{
! ~% b ` q' X1 e( _- y\\" K: f% m
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],
) \2 K: M; a% E1 Y3 G
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
7 o2 i5 b7 D, K' |( m; g; m6 z
};
. Y! ]; }6 i! I. W\\" T4 K& Z
[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))))); / @0 ?. r, A' c\" \0 h4 R, d
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:9 Q0 Z5 ~- B; r/ M7 X5 {
t=sys::clock();6 h& W/ R* Q% [5 p* [
s=0,x=0, . M5 K7 l0 [1 P# g
while{x<=1, //while循环算法； . \. y$ r3 \; L4 o' ~0 t% g1 V. L
y=1,
/ m, L6 j5 x) A\" [7 T( D
while{y<=2, ; n2 t* ~/ H: V7 g
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 7 I/ ?( f/ \7 p8 O0 q( g
y=y+0.0009 3 I, d5 c2 A5 n' E) p
},
7 u* U8 T; H. a* Y+ J1 V5 d\" ?
x=x+0.0011 ' `$ r) X# ?. z( C
}, ) F; J- ]7 @( o# J
s;' w! _; t. p. H. W
[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 y! o% r ?4 j7 w$ h' Zqbist 发表于 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"];$ g) }9 i- ^+ V7 g0 S2 k0 n7 m; }
(:a,b,k,t0)=+ Z5 u! u% U4 [0 O D4 Y
oo{
6 y+ I, _/ f. ] J \2 V7 H
a=rand[1000,1000], b=rand[1000,1000],
5 g3 @& D4 G/ I/ a$ K9 I
t0=clock(),
1 X' V% [5 m4 [
k=a*b, //矩阵乘
- w5 A4 U+ b: F+ Y' {! A/ E
k[1,3:5,9].outm()9 M& Z1 A0 A: g
},/ g( `3 o7 W\\" _ O/ X
[clock()-t0]/1000;3 ?. q( B, p* y# h* y

结果：

238.447 247.837 247.065 248.105 247.058
6 P) u. w, r& ~- r% d
244.123 249.925 247.553 243.981 250.016
, }# r. g2 |' {4 e
236.387 252.025 245.651 248.866 248.866% {( C+ G h& @( M9 s& O# L
2.219 秒' Q4 P. v4 [. l6 b

复制代码

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

!using["math","sys","XSLSF"];9 Z' s# X. P0 @$ m1 P' t8 E
(:a,b,k,t0)=& V8 d. ?+ [2 ?6 }% r) K# {8 M- x
oo{: k2 {; @, ~2 f2 W( |$ `: u6 c
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],& P( C8 t8 X0 X, D; V8 j7 Q
t0=clock(), T5 V( o* f2 V% ]0 H
rmul[k:a,b], //矩阵乘$ D6 i5 V3 b/ u, a3 S3 w) V
k[1,3:5,9].outm()1 s. l7 d5 s* R3 M0 O6 W
},$ r* o; Q0 c' ~; @+ Z) H3 s
[clock()-t0]/1000;
/ ^8 r% S% z7 z! e

结果：

262.121 247.583 260.529 259.548 258.328& n9 ^% `& p* s; \; s0 x
255.413 246.563 254.356 250.548 251.509! g1 C; @( f7 ^! Y) e5 t, \4 }' U
256.152 247.725 259.444 250.827 249.816
/ \$ |7 |( |1 j+ O
10.563 秒' v% \' N, [4 }+ \' @

复制代码

我也说一句

发表

使用道具举报

返回列表

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

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

浏览过的版块

新人进步奖

QQ

电话咨询

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