查看: 5451|回复: 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 |9 |; x' K( a9 Q a2 G; u
mvar:
\\" h, E5 e: s+ k0 ~9 F9 H( x* R1 C9 ?
oo{ //一般在oo函数中调用FcMath函数: i0 i- ]+ @- ~- s
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化
' F\\" a; _8 X, R
a.outm(), //输出矩阵a
# J\\" K1 |; }: u! S. ` C\\" W
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出
! V) T0 X1 G; z* w1 V: i
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出
3 T7 C7 E. }2 m9 r8 S7 |
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出
% V4 J; W- T: v2 C
}; c. Y) N, d) X5 j5 k

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.852615
, H# [0 Q5 B# b% i! d
0.630646 0.927048 0.440308 0.162857 0.556854: T. S* O$ R0 `
0.43309 0.34552 0.563919 0.937164 0.209641: o7 A6 w1 z( M/ R0 Z& `
0.603271 0.727676 0.130951 5.35736e-002 0.197937
8 D' g% r8 j1 J* \1 e4 a
0.576004 0.747589 1.17645e-002 0.363892 0.280777
6 P+ F, i7 ]1 h& x) L
0.646454 0.381088 0.58551 0.26387 0.93692
- z9 _4 ~) s, |6 @8 O) W
\! n E1 y0 S: X2 V& j- i
0.153259/ c: W4 I$ ~2 h% `9 u) k
0.162857
% q- w6 |) Q8 I6 X4 @6 Y\" V
0.937164
8 x+ H, W, u5 {/ M5 B; N' A& e
5.35736e-0026 i5 M& Q T+ J* D2 X- Q/ \
0.363892
2 [9 E: a1 N7 R+ p+ h
0.26387
# s) J4 X* U' C$ N2 _0 ?+ f) {$ {* _6 m
* J0 ~: l: T\" f- t; ]9 P) Q
0.603271 0.727676 0.130951 5.35736e-002 0.197937
u) V; R* `, G
; ? {; M, W8 Q9 u! U
0.130951 5.35736e-002
S; r- [6 K9 B+ I ?4 \
1.17645e-002 0.363892
. @2 K( F* K( }1 [6 G% J' c
0.58551 0.263871 J: @\" e. t1 G, i$ ?4 ~6 T
! ] j E9 d9 h) G1 ?* Y/ W

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义
0 X, d6 v\\" S+ D, N$ V; z0 h: x
{0 o$ R8 v! r+ }% v1 K
y1=x1*x1+x2*x2+x3*x3-1.0,+ d: p4 p7 @: e- U/ U: p
y2=2.0*x1*x1+x2*x2-4.0*x3,# Z\\" y% b$ Q6 U) t) F
y3=3.0*x1*x1-4.0*x2+x3*x3\\" n\\" [& M: s5 r1 C1 ?1 i
};( f: w: t+ W: Y: ^8 Q9 c. v5 t L
!using["math","sys"];9 ^9 g6 e/ r\\" S
mvar:
; V\\" l& i+ i! X) V5 I% t2 W; f
oo{
0 L! c: E! P ^! ^
x=array(3),
- \0 {( Q8 d) R\\" `, y
x.SA[0 : 1,1,1], //设置初值为1,1,1$ d: @$ \* A- g, k! r6 I( q+ K
i=netn[HFor("f"),x], //拟牛顿法解方程* H; ~, k5 I( A( o
x.outm(), //输出结果
1 n5 O& D( `3 |* S
i //返回迭代次数
2 A4 k% y- f! t
};! c1 Q\\" e/ b4 d3 z6 G

f(x1,x2,x3,y1,y2,y3)=      //函数定义 
0 X, d6 v\\" S+ D, N$ V; z0 h: x
{0 o$ R8 v! r+ }% v1 K

y1=x1*x1+x2*x2+x3*x3-1.0,+ d: p4 p7 @: e- U/ U: p

y2=2.0*x1*x1+x2*x2-4.0*x3,# Z\\" y% b$ Q6 U) t) F

y3=3.0*x1*x1-4.0*x2+x3*x3\\" n\\" [& M: s5 r1 C1 ?1 i

};( f: w: t+ W: Y: ^8 Q9 c. v5 t  L

!using["math","sys"];9 ^9 g6 e/ r\\" S

mvar: 
; V\\" l& i+ i! X) V5 I% t2 W; f
oo{ 
0 L! c: E! P  ^! ^
  x=array(3), 
- \0 {( Q8 d) R\\" `, y
  x.SA[0 : 1,1,1],       //设置初值为1,1,1$ d: @$ \* A- g, k! r6 I( q+ K

i=netn[HFor("f"),x],   //拟牛顿法解方程* H; ~, k5 I( A( o

x.outm(),              //输出结果 
1 n5 O& D( `3 |* S
  i                      //返回迭代次数 
2 A4 k% y- f! t
};! c1 Q\\" e/ b4 d3 z6 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$ O, [, h. G% [, j: h! u
clc
. O0 u2 s3 _- {5 \0 f
tic
k = zeros(5,5); % //生成5×5全0矩阵1 t( l2 r4 ^$ T
% 循环计算以下程序段100000次：3 r. F3 z! n; E4 F/ l* G% }; P# n5 `
for m = 1:1000003 c( H. O. K8 u5 j) A1 f3 `: B
a = rand(5,7);
2 v2 P7 G$ W5 ]) j1 G+ r
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
5 }\" d+ Q# {% D0 W9 }. x; ?
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
; f9 y+ C3 d; `. ^
end% Y1 m3 v$ W# `1 e
k
\" E& r# C6 ?+ o7 u2 D+ b
toc+ D* ?; a( J0 K) e5 q: n

复制代码

Forcal代码：

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

!using["math","sys"];
9 X! K\\" ?5 o! L\\" ] a
mvar:7 u\\" Q! H9 V' d \
t0=clock(),
# a' j M9 g2 Y* M
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0: D. F$ ~8 c' l( G7 W
i=0,(i<1000 00).while{ //循环计算1000 00次
o7 F; G6 H9 A4 q$ j9 u
oo{, B1 I\\" f) q' a+ t, T
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
$ \8 c0 A! r8 S% t4 x, ^
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)/ @+ J/ ~5 H+ s5 k! s) C6 s
},; W' T/ d8 s/ k7 J0 B\\" G+ J1 }
i++8 ?# B' c* D4 z- j1 V
},
: N# J$ B6 _: {+ \3 }' b/ g
k.outm(), //输出矩阵k，然后销毁k/ L& F+ J2 s) w% O6 y
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"]; 
9 X! K\\" ?5 o! L\\" ]  a
mvar:7 u\\" Q! H9 V' d  \

t0=clock(), 
# a' j  M9 g2 Y* M
oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0: D. F$ ~8 c' l( G7 W

i=0,(i<1000 00).while{ //循环计算1000 00次 
  o7 F; G6 H9 A4 q$ j9 u
  oo{, B1 I\\" f) q' a+ t, T

a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化 
$ \8 c0 A! r8 S% t4 x, ^
    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)/ @+ J/ ~5 H+ s5 k! s) C6 s

},; W' T/ d8 s/ k7 J0 B\\" G+ J1 }

i++8 ?# B' c* D4 z- j1 V

}, 
: N# J$ B6 _: {+ \3 }' b/ g
k.outm(),             //输出矩阵k，然后销毁k/ L& F+ J2 s) w% O6 y

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

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

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

!using["math","sys"];
+ Q7 \1 |5 Y8 c& o9 `$ H8 C! N
(:t0,k,i,a,b)=! m( ~3 i% ?' P- \; Z6 |
{, W. h# ~: ]* o9 F7 g1 W( w4 Y
t0=clock(),
2 W$ Q3 V4 }0 L4 c0 s
k=zeros[5,5],
/ S9 |\\" H\\" ^+ v* r* Y+ I
i=0,(i<1000 00).while{7 E* F: A1 k7 m# u' h8 b
oo{
3 D0 u! G7 f7 D/ h' f+ Q# ]
a=rand[5,7], b=rand[7,5],
6 B! }* x* v5 l- Y! Z' C
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
, r' T+ z2 K7 {, ?3 ~
},9 A0 u* m/ s1 s3 }: R% f$ y
i++
3 L/ ^7 k6 v2 @, w+ K
},0 k* O% d% }6 q0 L, z\\" r; ^% t* C
k.outm().delete(),: A& N* Z1 Z' g) a
[clock()-t0]/1000
0 T ^: m# b/ Z4 N5 M+ I
};

!using["math","sys"]; 
+ Q7 \1 |5 Y8 c& o9 `$ H8 C! N
(:t0,k,i,a,b)=! m( ~3 i% ?' P- \; Z6 |

{, W. h# ~: ]* o9 F7 g1 W( w4 Y

t0=clock(), 
2 W$ Q3 V4 }0 L4 c0 s
  k=zeros[5,5], 
/ S9 |\\" H\\" ^+ v* r* Y+ I
  i=0,(i<1000 00).while{7 E* F: A1 k7 m# u' h8 b

oo{ 
3 D0 u! G7 f7 D/ h' f+ Q# ]
      a=rand[5,7], b=rand[7,5], 
6 B! }* x* v5 l- Y! Z' C
      k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
, r' T+ z2 K7 {, ?3 ~
    },9 A0 u* m/ s1 s3 }: R% f$ y

i++ 
3 L/ ^7 k6 v2 @, w+ K
  },0 k* O% d% }6 q0 L, z\\" r; ^% t* C

k.outm().delete(),: A& N* Z1 Z' g) a

[clock()-t0]/1000 
0 T  ^: m# b/ Z4 N5 M+ I
};

在我的电脑上运行时间为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"];
. r1 e$ K3 K+ t+ Z6 F A' N! D% q
mvar:
* Q$ d+ v9 q1 B
oo{
9 c) Q& I, r3 O) k0 o- J! h
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],0 {3 ^* F- u3 O7 v: Z5 u2 S. T
a1.outm[5,1,1],: l- _4 V5 B$ D6 }+ P
a2.outm[5,1,1],0 o4 ~+ ]8 y0 E6 \* O0 H& ]8 g% a
a3.outm[5,1,1],! K$ p1 l\\" R: L1 H. l
a4.outm[5,1,1],
/ d) Q( P. J. P
a=a1+a2+a3+a4,( |8 S& E# J4 r\\" |0 V+ s$ m
a.outm[5,1,1],
S: d) H3 p( |2 o8 j, B: w
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]
( d) u5 e& ?' y3 y O. n7 i! m\\" h
};
; W( [7 T; I! O, l/ b; |\\" _

!using["math"]; 
. r1 e$ K3 K+ t+ Z6 F  A' N! D% q
mvar: 
* Q$ d+ v9 q1 B
oo{ 
9 c) Q& I, r3 O) k0 o- J! h
  ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],0 {3 ^* F- u3 O7 v: Z5 u2 S. T

a1.outm[5,1,1],: l- _4 V5 B$ D6 }+ P

a2.outm[5,1,1],0 o4 ~+ ]8 y0 E6 \* O0 H& ]8 g% a

a3.outm[5,1,1],! K$ p1 l\\" R: L1 H. l

a4.outm[5,1,1], 
/ d) Q( P. J. P
  a=a1+a2+a3+a4,( |8 S& E# J4 r\\" |0 V+ s$ m

a.outm[5,1,1], 
  S: d) H3 p( |2 o8 j, B: w
  Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[] 
( d) u5 e& ?' y3 y  O. n7 i! m\\" h
}; 
; W( [7 T; I! O, l/ b; |\\" _

说明：
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"];
' X f\\" `# a& Q\\" _; U
mvar:/ c8 Q. R2 q: f5 ~
(:p1,p2,p3,a,b)=, |9 x# _$ b5 U0 t$ p. Y
{* u9 l\\" p7 `6 m1 n
oo{7 y+ k+ x, r- [5 B6 _+ j$ m7 i! |
a=array[1000].rand(),0 s' @+ ?; e% P* \( v9 _( y
b=array[1000].rand(),% A' W) n/ p! q1 q; X7 r3 D ]
p1=array[1000,1000],
1 D ~; h/ o8 X+ v( D7 E
p2=array[1000,1000],( d, S7 ~0 r1 n9 a. F. @5 i
p3=array[1000,1000],
8 D* J; K' G% P) m3 z, i
t0=clock(),$ `( \( d0 |5 V% F/ t, W1 U2 B
ndgrid(a,b,&A,&B),6 t3 a9 J: p! X' w+ {) R
p1.=A+B
' P. j! s, J! t' e
},
! q2 ~+ r# B, T0 n& Q\\" y, P( L
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},5 b. P5 I! h* g
lena=FCDLen(a),( L( c/ z' M. \+ i' }3 W
lenb=FCDLen(b),
$ } b$ { ^' S% N+ l4 q- {2 |
t0=clock(),4 L/ ]1 g2 u& X# I+ O* t( N. x& v
m = lenb-1, (m>=0).while{
* C4 D0 c( Z\\" ~# v
oo{p2(m,neg) = a+rn[b(m)]},' B8 A! s8 g) L2 ?( w) ]
m--9 i' i9 O( G/ M. t: A
},
3 ]$ `* n7 g1 ~# R\\" |$ D
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},. P9 l\\" `5 B0 q9 c& b
t0=clock(),5 p# s5 o- ~* h* n3 v; s
m = lenb-1, (m>=0).while{( L3 ]0 s/ I5 V2 b
n = lena-1, (n>=0).while{
$ e1 L8 H8 D; o
//p3(m,n) = a(n)+b(m), //用这句还要慢一些 T& H+ \ z& c. `9 g# P, A& w6 R
A(p3,m,n) = A(a,n)+A(b,m),
\\" J l& s& C( ~% \7 Q. J. Y4 M( o
n--3 F. D2 x\\" E8 D# F$ k( M1 I
},
* J8 K: e1 q+ a$ I S2 i
m--6 k) J$ l( A: b. j$ Q
},. ?/ |. }! F4 E1 h. Q2 i* ~ e
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}
( A v/ Z/ Z$ Y' t' ^, o0 C, t
};
4 ~+ R) {, t( d6 F3 k6 A

!using["math","sys"]; 
' X  f\\" `# a& Q\\" _; U
mvar:/ c8 Q. R2 q: f5 ~

(:p1,p2,p3,a,b)=, |9 x# _$ b5 U0 t$ p. Y

{* u9 l\\" p7 `6 m1 n

oo{7 y+ k+ x, r- [5 B6 _+ j$ m7 i! |

a=array[1000].rand(),0 s' @+ ?; e% P* \( v9 _( y

b=array[1000].rand(),% A' W) n/ p! q1 q; X7 r3 D  ]

p1=array[1000,1000], 
1 D  ~; h/ o8 X+ v( D7 E
    p2=array[1000,1000],( d, S7 ~0 r1 n9 a. F. @5 i

p3=array[1000,1000], 
8 D* J; K' G% P) m3 z, i
    t0=clock(),$ `( \( d0 |5 V% F/ t, W1 U2 B

ndgrid(a,b,&A,&B),6 t3 a9 J: p! X' w+ {) R

p1.=A+B 
' P. j! s, J! t' e
  }, 
! q2 ~+ r# B, T0 n& Q\\" y, P( L
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},5 b. P5 I! h* g

lena=FCDLen(a),( L( c/ z' M. \+ i' }3 W

lenb=FCDLen(b), 
$ }  b$ {  ^' S% N+ l4 q- {2 |
  t0=clock(),4 L/ ]1 g2 u& X# I+ O* t( N. x& v

m = lenb-1, (m>=0).while{ 
* C4 D0 c( Z\\" ~# v
    oo{p2(m,neg) = a+rn[b(m)]},' B8 A! s8 g) L2 ?( w) ]

m--9 i' i9 O( G/ M. t: A

}, 
3 ]$ `* n7 g1 ~# R\\" |$ D
  printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},. P9 l\\" `5 B0 q9 c& b

t0=clock(),5 p# s5 o- ~* h* n3 v; s

m = lenb-1, (m>=0).while{( L3 ]0 s/ I5 V2 b

n = lena-1, (n>=0).while{ 
$ e1 L8 H8 D; o
      //p3(m,n) = a(n)+b(m),  //用这句还要慢一些  T& H+ \  z& c. `9 g# P, A& w6 R

A(p3,m,n) = A(a,n)+A(b,m), 
\\" J  l& s& C( ~% \7 Q. J. Y4 M( o
      n--3 F. D2 x\\" E8 D# F$ k( M1 I

}, 
* J8 K: e1 q+ a$ I  S2 i
    m--6 k) J$ l( A: b. j$ Q

},. ?/ |. }! F4 E1 h. Q2 i* ~  e

printff{"\r\nfor2: {1,r}",[clock()-t0]/1000} 
( A  v/ Z/ Z$ Y' t' ^, o0 C, t
}; 
4 ~+ R) {, t( d6 F3 k6 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"]; n4 D7 d% c+ P: q% S6 k4 K* ?
mvar:) `# B4 ~! H' m5 _8 O& f8 H5 |
t=clock(),
- s! R0 Q H0 [* ^! F
oo{4 F) }) k( N, v* w0 `4 g+ l
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],- V, B# R( s4 t6 |' s5 [: T
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]- ?8 L* E\\" ]! c! P* k; h
};
6 W3 B9 X/ s; x8 S
[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 ]/ p\" ?5 I7 |* _, q. Z
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:8 [8 O% P0 ?9 a$ g7 D
t=sys::clock();6 b6 t5 }8 ], M\" P
s=0,x=0, 6 F6 u) y7 s. q A, }8 p
while{x<=1, //while循环算法；
0 d$ t* u6 i; D
y=1, , E: g5 `$ m3 ^
while{y<=2,
. T* v* Z! ^ Q
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))), 3 S6 i8 ]& f2 n* o& D M8 _$ v7 y
y=y+0.0009 ) ]% O. y7 t2 G/ _- E, l
}, Q\" T: `4 Z+ p. R
x=x+0.0011
+ @( W: g8 W: l
}, ) s6 q4 k, e. D4 W3 S* ?
s;
/ \- a! [3 R0 w4 k1 _- G& j
[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 编辑

好深奥！~~~~3 [" Q$ \7 E& j7 T7 h3 l0 U
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"];6 G9 i7 X1 j; F& ?6 N
(:a,b,k,t0)=
8 M4 A @5 a( m8 a% I1 A\\" i
oo{
% I# k7 n$ K7 x- L- y- R. Z. O' [
a=rand[1000,1000], b=rand[1000,1000],8 J+ @: n; g3 h3 i. }. M
t0=clock(),\\" T! Q& f1 p% w* ^# x6 D
k=a*b, //矩阵乘
( A# s' ?/ F+ G* \9 T
k[1,3:5,9].outm()8 c/ Z( ~# s% {4 i* u! ~
},* Z* a, }/ e& o9 t+ x
[clock()-t0]/1000;$ }4 F) H# o# [- I

结果：

238.447 247.837 247.065 248.105 247.0585 r4 z: a0 m& c1 h0 Q* T/ _9 {
244.123 249.925 247.553 243.981 250.016
9 F/ o5 C6 p$ h* t\" H; X\" H6 F1 w
236.387 252.025 245.651 248.866 248.866! v3 K: f- f3 O
2.219 秒
. P! L1 {. A5 A0 t% @6 ~* Q

复制代码

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

!using["math","sys","XSLSF"];$ V$ ~& W% C9 j# V# J! n
(:a,b,k,t0)=\\" S' T\\" p# O0 j0 Q
oo{
% y3 P* ]2 k; |2 r7 U j, l
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],0 y/ g9 Q9 y( @- | B( ?$ p
t0=clock(),3 { [ f1 n2 D/ N& W0 Q6 m
rmul[k:a,b], //矩阵乘
* t P# L' J# k\\" U( ~
k[1,3:5,9].outm(): z/ ^. N- M9 s. o ?1 G& i
},; z( N0 F8 @9 j\\" j7 n, r, I
[clock()-t0]/1000;
* O' R$ S& X' R0 B) s) O

结果：

262.121 247.583 260.529 259.548 258.328
! K- d/ i) b2 ?. V; C4 g, C6 O7 V+ q
255.413 246.563 254.356 250.548 251.5091 u6 j0 h6 ^: F7 R
256.152 247.725 259.444 250.827 249.8167 `) R0 b5 z' d/ _$ Y
10.563 秒1 M# p) M3 o! X7 [* S1 E; ]: ?

复制代码

我也说一句

发表

使用道具举报

返回列表

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

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

新人进步奖

QQ

电话咨询

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