查看: 5732|回复: 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 m& ~ u6 _\\" |) i
mvar:\\" r( I4 c8 H L% Q/ h
oo{ //一般在oo函数中调用FcMath函数
, y: o8 L# c+ Z) Z( f0 b6 Z
a=rand[6,5], //生成6×5矩阵a，用0～1之间随机数初始化
# T5 [ S/ G3 S) J7 E
a.outm(), //输出矩阵a
$ o3 G( H$ k$ {. y% `
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出, V u, X$ b1 h2 [
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出' e7 M2 T& ~( ^8 d' n) J6 t! F& c
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出
4 a8 J* h( m' P) l. F
};
6 O% K6 J F! U* M2 F& y/ I- U

结果：

0.211319 4.91638e-002 0.144638 0.153259 0.852615' A) M1 m/ ]( a4 X6 ]; }
0.630646 0.927048 0.440308 0.162857 0.556854- i) Q) ^% H% A' J# S) F3 v, X
0.43309 0.34552 0.563919 0.937164 0.209641
4 v\" c) a$ P7 Z( u$ t8 u4 a
0.603271 0.727676 0.130951 5.35736e-002 0.197937
6 m, @: O$ }3 s: T( a3 i! [
0.576004 0.747589 1.17645e-002 0.363892 0.280777- C, c1 @; j( O- b0 f/ F6 M: q
0.646454 0.381088 0.58551 0.26387 0.93692
% c: ]5 ^, Z9 w$ s- `4 @
5 d* Y8 ?) @% L8 o2 N: p
0.153259
3 F) x) b3 ^' d) f+ S
0.1628571 r' m% V4 o2 H
0.937164
; H8 h/ }. q) E! \0 |/ k
5.35736e-002
/ F& P0 r0 ?9 S' ^- {: i+ K
0.3638920 T- \) H) W1 I1 ^+ Y5 v9 z. @
0.26387# O' m$ s& L) q1 d$ N/ N& W
/ r1 a\" K( F( Y+ W
0.603271 0.727676 0.130951 5.35736e-002 0.1979373 P* z0 t* o! u5 Y5 b: n) u
) @6 Y }& @: |% }3 T; A. j
0.130951 5.35736e-002) Z; F3 d. C. ~% x4 c9 `, L3 q1 |( F
1.17645e-002 0.363892
. v( u7 Q- a2 t( o8 i# E' ?8 V1 L: e
0.58551 0.26387
, Q9 g. `2 U5 c2 _4 \1 o# T
: k% W9 e8 n8 T/ t7 d, ~1 d

复制代码

例子2代码：

f(x1,x2,x3,y1,y2,y3)= //函数定义0 E m6 O- Q. x, |# ^+ Y+ b
{
, O8 t2 P: \4 y: ?
y1=x1*x1+x2*x2+x3*x3-1.0,, n0 L# f, O m, Z
y2=2.0*x1*x1+x2*x2-4.0*x3,
$ Z5 ]$ O1 ]' w) P1 P; d
y3=3.0*x1*x1-4.0*x2+x3*x3
, T: ~8 o9 y( ]3 c# f1 x
};9 D/ i, |5 f\\" g
!using["math","sys"];/ \ e3 Y5 ?; J& B
mvar:1 }4 i+ A& g7 l\\" {; c y- _
oo{
% L, `3 U. b0 n
x=array(3),
# @: j5 |1 D+ L+ m
x.SA[0 : 1,1,1], //设置初值为1,1,1; T0 C2 K+ N+ n9 B& h+ r
i=netn[HFor("f"),x], //拟牛顿法解方程
- l4 e+ [! r; o( M
x.outm(), //输出结果. A4 B. s\\" v; i1 E
i //返回迭代次数
4 T9 a5 n# R* q7 t7 |
};. s! c* {+ t$ o# E- b3 X

f(x1,x2,x3,y1,y2,y3)=      //函数定义0 E  m6 O- Q. x, |# ^+ Y+ b

{ 
, O8 t2 P: \4 y: ?
    y1=x1*x1+x2*x2+x3*x3-1.0,, n0 L# f, O  m, Z

y2=2.0*x1*x1+x2*x2-4.0*x3, 
$ Z5 ]$ O1 ]' w) P1 P; d
    y3=3.0*x1*x1-4.0*x2+x3*x3 
, T: ~8 o9 y( ]3 c# f1 x
};9 D/ i, |5 f\\" g

!using["math","sys"];/ \  e3 Y5 ?; J& B

mvar:1 }4 i+ A& g7 l\\" {; c  y- _

oo{ 
% L, `3 U. b0 n
  x=array(3), 
# @: j5 |1 D+ L+ m
  x.SA[0 : 1,1,1],       //设置初值为1,1,1; T0 C2 K+ N+ n9 B& h+ r

i=netn[HFor("f"),x],   //拟牛顿法解方程 
- l4 e+ [! r; o( M
  x.outm(),              //输出结果. A4 B. s\\" v; i1 E

i                      //返回迭代次数 
4 T9 a5 n# R* q7 t7 |
};. s! c* {+ t$ o# E- b3 X

结果：
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
' x# k0 @1 z* I1 V
clc
& q+ O2 p2 U- `
tic3 [& n. _8 S1 X1 f
k = zeros(5,5); % //生成5×5全0矩阵
+ Y) d3 B) M& Q& u
% 循环计算以下程序段100000次：. K3 |, o2 E! ]; w. B' u! f$ v; P
for m = 1:100000) S% r7 @9 Z/ G! S
a = rand(5,7);
9 w5 s' x+ {5 H& w6 x$ G
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化9 @7 s: G9 |# `; g; T7 O& I6 ~
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
: U# d) L7 b& V\" O1 R( i
end) X H2 p; {! B! j }) f
k
% B, R$ E\" Y Q5 x; \
toc% l9 x5 Y3 y\" L. F/ Z. |

复制代码

Forcal代码：

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

!using["math","sys"];( T) @9 u! b n7 u. C7 q
mvar:
. s7 O2 M6 x# T6 U. A
t0=clock(),4 h# [, M. J\\" u1 g\\" O
oo{k=zeros[5,5]}, //生成5×5矩阵k，初始化为0( h- g) x! z$ E( r
i=0,(i<1000 00).while{ //循环计算1000 00次, E0 r8 [+ T4 f1 g9 g9 l
oo{) R/ n, m# Z0 T2 M$ V
a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化9 I i4 f0 [0 n! v7 K' g
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)& z. W+ ~! b* W* V
},
$ d9 Y, Y: f4 j y4 r: g
i++
; W9 x! f' g5 J% m
},' `* O; E! @1 v; R& [
k.outm(), //输出矩阵k，然后销毁k0 r% A& z F+ K, } y9 _. y0 ?: q* n
[clock()-t0]/1000; //得到计算时间，秒

!using["math","sys"];( T) @9 u! b  n7 u. C7 q

mvar: 
. s7 O2 M6 x# T6 U. A
t0=clock(),4 h# [, M. J\\" u1 g\\" O

oo{k=zeros[5,5]},     //生成5×5矩阵k，初始化为0( h- g) x! z$ E( r

i=0,(i<1000 00).while{ //循环计算1000 00次, E0 r8 [+ T4 f1 g9 g9 l

oo{) R/ n, m# Z0 T2 M$ V

a=rand[5,7], b=rand[7,5], //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化9 I  i4 f0 [0 n! v7 K' g

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)& z. W+ ~! b* W* V

}, 
$ d9 Y, Y: f4 j  y4 r: g
  i++ 
; W9 x! f' g5 J% m
},' `* O; E! @1 v; R& [

k.outm(),             //输出矩阵k，然后销毁k0 r% A& z  F+ K, }  y9 _. y0 ?: q* n

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

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

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

!using["math","sys"]; u j9 ~ N9 [2 \8 K' h
(:t0,k,i,a,b)=5 c2 N3 t' A& e: t& [\\" o
{, R- C* [) |+ s, J1 ?' ~
t0=clock(),
! ^: r+ {% s, [) B4 J) I, x5 W- R
k=zeros[5,5],5 E' t! i5 ^* ?+ ^- n
i=0,(i<1000 00).while{
, x7 x0 T* U4 ^# j) m6 l# v
oo{6 S& x1 A$ [; e9 {: F3 f+ @# T
a=rand[5,7], b=rand[7,5],$ l, H9 E! v7 A+ Z3 @
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
9 \& Y: y, M2 r9 H' K
},
6 }4 h\\" g2 p( }2 q6 a5 z
i+++ ` p! z5 g$ q1 w
},\\" q& G$ ~; Q( }+ J\\" X\\" p# h
k.outm().delete(),$ q8 M4 X# s+ h- J7 \0 [. f
[clock()-t0]/1000- A' E( u, a4 y% H0 A% J+ l
};

!using["math","sys"];  u  j9 ~  N9 [2 \8 K' h

(:t0,k,i,a,b)=5 c2 N3 t' A& e: t& [\\" o

{, R- C* [) |+ s, J1 ?' ~

t0=clock(), 
! ^: r+ {% s, [) B4 J) I, x5 W- R
  k=zeros[5,5],5 E' t! i5 ^* ?+ ^- n

i=0,(i<1000 00).while{ 
, x7 x0 T* U4 ^# j) m6 l# v
    oo{6 S& x1 A$ [; e9 {: F3 f+ @# T

a=rand[5,7], b=rand[7,5],$ l, H9 E! v7 A+ Z3 @

k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg) 
9 \& Y: y, M2 r9 H' K
    }, 
6 }4 h\\" g2 p( }2 q6 a5 z
    i+++ `  p! z5 g$ q1 w

},\\" q& G$ ~; Q( }+ J\\" X\\" p# h

k.outm().delete(),$ q8 M4 X# s+ h- J7 \0 [. f

[clock()-t0]/1000- A' E( u, a4 y% H0 A% J+ l

};

在我的电脑上运行时间为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"];
. q( t2 r$ K1 w
mvar:) x( Z6 d7 p' e& o! O. d
oo{
& L$ o) s F\\" X7 L6 d, y0 h+ P
ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],7 H- L, S# _/ p( Z4 n% z2 e* S. t
a1.outm[5,1,1],
+ B8 c3 p; S' r& }/ @& G( Y; |1 `
a2.outm[5,1,1],
) U% T3 _7 }\\" f5 u0 w0 o
a3.outm[5,1,1],
Z) v- v# S* }1 ]% y
a4.outm[5,1,1],; x' X* Z( d: L4 n* q8 {$ V
a=a1+a2+a3+a4,
& @* G$ Z9 x9 `# [, S* V* E- d
a.outm[5,1,1],/ V, V\\" S& c) Q6 N3 a! U- h0 s$ j
Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]7 T. y5 ]4 }( }3 O
};
: v& i5 \/ K+ H6 a- k

!using["math"]; 
. q( t2 r$ K1 w
mvar:) x( Z6 d7 p' e& o! O. d

oo{ 
& L$ o) s  F\\" X7 L6 d, y0 h+ P
  ndgrid[linspace(1,2,2),linspace(3,5,3),linspace(6,7,2),linspace(8,12,5),&a1,&a2,&a3,&a4],7 H- L, S# _/ p( Z4 n% z2 e* S. t

a1.outm[5,1,1], 
+ B8 c3 p; S' r& }/ @& G( Y; |1 `
  a2.outm[5,1,1], 
) U% T3 _7 }\\" f5 u0 w0 o
  a3.outm[5,1,1], 
  Z) v- v# S* }1 ]% y
  a4.outm[5,1,1],; x' X* Z( d: L4 n* q8 {$ V

a=a1+a2+a3+a4, 
& @* G$ Z9 x9 `# [, S* V* E- d
  a.outm[5,1,1],/ V, V\\" S& c) Q6 N3 a! U- h0 s$ j

Sum[a].outm[5,1,1].Sum[].outm[5,1,1].Sum[].outm[5,1,1].Sum[]7 T. y5 ]4 }( }3 O

}; 
: v& i5 \/ K+ H6 a- k

说明：
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 k& Z6 D+ o& h! ?
mvar:
4 w P% P, r# g( H/ k4 O
(:p1,p2,p3,a,b)=$ f0 i9 g4 k ~7 Z, C5 n' U5 t: R
{3 U( f1 s# G( S; S% q+ ?
oo{
Y5 A9 Y) P7 d0 {( r
a=array[1000].rand(),& x/ x5 \/ E2 Q' Z1 {4 J8 z; M2 h
b=array[1000].rand(),# [! n, R4 [9 z: H8 ?
p1=array[1000,1000],
7 V$ b* B4 \& ~% ?9 o: _) H% x* u
p2=array[1000,1000],
8 F6 T$ X\\" r( j* T
p3=array[1000,1000],; O3 [- G. j' o9 X1 w) K- Q9 S$ \
t0=clock(),8 S7 s5 v3 d- d4 b
ndgrid(a,b,&A,&B),
5 P0 E3 p( P. K( D2 z O$ [8 a
p1.=A+B
! c0 o$ }9 G( ^* S\\" G2 W
},
4 |/ y2 }$ j1 m6 s; A( @' @( |4 h
printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},+ r1 T# K\\" D8 Y9 _* G
lena=FCDLen(a),
7 u6 ?9 s) _2 P) V* `
lenb=FCDLen(b),
# \* w5 A5 Z2 c' l6 d5 s; O
t0=clock(),
# \6 X. r' a2 ?0 K& }( ~/ P
m = lenb-1, (m>=0).while{
0 S! H3 e3 f\\" H* `, Y$ q
oo{p2(m,neg) = a+rn[b(m)]},+ r5 \4 L- _, ?) y& B
m--3 e2 U# i+ R5 g/ q: R
},2 u$ f& k+ W. P6 s* o5 M1 n
printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},\\" b* M: _- P4 \! a
t0=clock(),! e- v6 K4 ^, |' J/ ~- c/ ]: P; G+ }
m = lenb-1, (m>=0).while{* p# A! X1 k* Y1 P9 N. P
n = lena-1, (n>=0).while{
% v* w$ {: I* `% {! X' p7 R
//p3(m,n) = a(n)+b(m), //用这句还要慢一些; {/ y+ ?6 I, J
A(p3,m,n) = A(a,n)+A(b,m),
9 k/ ]\\" J. q2 D1 V; I- B4 F4 h2 R
n--7 [' N+ B6 |* n! E- ~7 v k% K
},
H0 n6 a- g7 _0 n\\" z% R
m--2 k, N& S\\" K. d8 M- F# K
},
* ?8 ^* [2 M/ H; I9 ~+ G
printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}( k; ~4 Z: E( G9 L( C. k! L
};2 n$ b9 V* H) m6 A' e& s. ^

!using["math","sys"];5 k& Z6 D+ o& h! ?

mvar: 
4 w  P% P, r# g( H/ k4 O
(:p1,p2,p3,a,b)=$ f0 i9 g4 k  ~7 Z, C5 n' U5 t: R

{3 U( f1 s# G( S; S% q+ ?

oo{ 
  Y5 A9 Y) P7 d0 {( r
    a=array[1000].rand(),& x/ x5 \/ E2 Q' Z1 {4 J8 z; M2 h

b=array[1000].rand(),# [! n, R4 [9 z: H8 ?

p1=array[1000,1000], 
7 V$ b* B4 \& ~% ?9 o: _) H% x* u
    p2=array[1000,1000], 
8 F6 T$ X\\" r( j* T
    p3=array[1000,1000],; O3 [- G. j' o9 X1 w) K- Q9 S$ \

t0=clock(),8 S7 s5 v3 d- d4 b

ndgrid(a,b,&A,&B), 
5 P0 E3 p( P. K( D2 z  O$ [8 a
    p1.=A+B 
! c0 o$ }9 G( ^* S\\" G2 W
  }, 
4 |/ y2 }$ j1 m6 s; A( @' @( |4 h
  printff{"\r\nndgrid: {1,r}",[clock()-t0]/1000},+ r1 T# K\\" D8 Y9 _* G

lena=FCDLen(a), 
7 u6 ?9 s) _2 P) V* `
  lenb=FCDLen(b), 
# \* w5 A5 Z2 c' l6 d5 s; O
  t0=clock(), 
# \6 X. r' a2 ?0 K& }( ~/ P
  m = lenb-1, (m>=0).while{ 
0 S! H3 e3 f\\" H* `, Y$ q
    oo{p2(m,neg) = a+rn[b(m)]},+ r5 \4 L- _, ?) y& B

m--3 e2 U# i+ R5 g/ q: R

},2 u$ f& k+ W. P6 s* o5 M1 n

printff{"\r\nfor1: {1,r}",[clock()-t0]/1000},\\" b* M: _- P4 \! a

t0=clock(),! e- v6 K4 ^, |' J/ ~- c/ ]: P; G+ }

m = lenb-1, (m>=0).while{* p# A! X1 k* Y1 P9 N. P

n = lena-1, (n>=0).while{ 
% v* w$ {: I* `% {! X' p7 R
      //p3(m,n) = a(n)+b(m),  //用这句还要慢一些; {/ y+ ?6 I, J

A(p3,m,n) = A(a,n)+A(b,m), 
9 k/ ]\\" J. q2 D1 V; I- B4 F4 h2 R
      n--7 [' N+ B6 |* n! E- ~7 v  k% K

}, 
  H0 n6 a- g7 _0 n\\" z% R
    m--2 k, N& S\\" K. d8 M- F# K

}, 
* ?8 ^* [2 M/ H; I9 ~+ G
  printff{"\r\nfor2: {1,r}",[clock()-t0]/1000}( k; ~4 Z: E( G9 L( C. k! L

};2 n$ b9 V* H) m6 A' e& s. ^

结果：
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"];
6 }. A* R+ C# ?
mvar:8 o: K/ F& l! Z1 B; B
t=clock(),
' k3 i. a/ [- c- ~: [( N5 i
oo{+ w; r3 r6 i1 e5 D {( x
ndgrid[linspacex(0,1,0.0011),linspacex(1,2,0.0009),&x,&y],. _; A% q* M\\" }1 ^' H
Sum[Cos(rn(1)-Sin(rn(1.2)*x^(y/rn(2))+Cos(rn(1)-Sin(rn(1.2)*y^(x/rn (2)))))),0]
+ x9 c% r0 O& K, K. P0 o+ q
};
% g! K' b+ y3 F0 H
[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))))); ; |! g+ @; t( ?7 m1 z
sum["f",0,1,0.0011,1,2,0.0009];

复制代码

结果：
1008606.64947441
0.719 //时间

3、while循环

mvar:
/ @$ M% x, J5 z
t=sys::clock();2 H; G1 h; B' `\" X7 u7 S1 w/ o
s=0,x=0,
6 X0 t2 n S& R# m) O$ N
while{x<=1, //while循环算法； - [8 c3 e4 A) D\" H6 X: c8 ?' K3 C3 p* Z
y=1, I E; G- }5 P) Z
while{y<=2,
. r+ @5 Q4 `& P
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
6 L% R# I4 S9 r# c3 R: c; v2 \
y=y+0.0009
3 }$ B7 s6 S ?0 X2 Y+ X
},
/ l. |/ z1 ~8 j* d K$ T$ f
x=x+0.0011
$ x7 J& E# q0 j7 P/ n' R# `3 ^5 {, U
},
& V7 Z- ?0 U& r. {
s;
' c. s4 [! y2 w! O8 D
[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 编辑

好深奥！~~~~; a/ i+ N1 y* E1 n/ i& b
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"];
G: _; ^( K, F8 b4 u! q2 f. X
(:a,b,k,t0)=
) t5 n! u( E& [) |/ t' b9 v
oo{
' U- j2 f! c. J/ B\\" y+ T
a=rand[1000,1000], b=rand[1000,1000],
! b\\" i+ u) ?# n( j1 O1 \4 x2 \
t0=clock(),
Z9 c7 b1 @& A4 `
k=a*b, //矩阵乘
$ @( W& m! w9 j
k[1,3:5,9].outm()* m9 P3 I! W, K: P% v
},
& | m+ Y `% U% v
[clock()-t0]/1000;' i! S; I) m( V2 i b' q

结果：

238.447 247.837 247.065 248.105 247.058
\" ^% I, a, Q\" P q$ d) X
244.123 249.925 247.553 243.981 250.016
/ q$ n ^, d4 a4 v
236.387 252.025 245.651 248.866 248.866# g3 [/ [ |# C6 s5 O0 N* m
2.219 秒5 U5 F6 v, G9 V D( A- c

复制代码

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

!using["math","sys","XSLSF"];( T; i# U0 L, o! Y2 L
(:a,b,k,t0)=1 S7 B$ |5 A0 ]# ]$ n% T6 o
oo{
# m! m' ^ {* z6 [
a=rand[1000,1000], b=rand[1000,1000], k=array[1000,1000],3 y& J6 p1 ~! r, d% u$ }) ?
t0=clock(),# z4 u0 S: v. h
rmul[k:a,b], //矩阵乘
3 P* i\\" \& J! |/ O4 r* h
k[1,3:5,9].outm()5 I. K5 T! l3 X' Q, I$ {8 L( ` \
},; Y+ z2 J/ _8 V+ x7 J. S\\" X1 \' q
[clock()-t0]/1000;
0 m. Y3 O6 a! _

结果：

262.121 247.583 260.529 259.548 258.328
9 o( Q! b/ x E6 V5 N6 T8 b8 Y
255.413 246.563 254.356 250.548 251.509
) C4 [! Y: r3 T) m# q% s' M
256.152 247.725 259.444 250.827 249.816
5 z; F3 H& K8 t6 ~) [
10.563 秒
7 L: s: }4 t2 x; _( e# Q& G& w6 v

复制代码

我也说一句

发表

使用道具举报

返回列表

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

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

新人进步奖

QQ

电话咨询

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