Forcal数学库FcMath：以矩阵运算为基础
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 例子1代码：
!using["math"];
mvar:
oo{ //一般在oo函数中调用FcMath函数
a=rand, //生成6×5矩阵a，用0～1之间随机数初始化
a.outm(), //输出矩阵a
a.subg(neg:3).outm(), //取矩阵a第4列所有元素组成子矩阵，并输出
a.subg(3:neg).outm(), //取矩阵a第4行所有元素组成子矩阵，并输出
a.subg(3,5:2,3).outm() //取矩阵a第4～6行，3～4列所有元素组成子矩阵，并输出
};
结果：
0.211319 4.91638e-002 0.144638 0.153259 0.852615
0.630646 0.927048 0.440308 0.162857 0.556854
0.43309 0.34552 0.563919 0.937164 0.209641
0.603271 0.727676 0.130951 5.35736e-002 0.197937
0.576004 0.747589 1.17645e-002 0.363892 0.280777
0.646454 0.381088 0.58551 0.26387 0.93692
0.153259
0.162857
0.937164
5.35736e-002
0.363892
0.26387
0.603271 0.727676 0.130951 5.35736e-002 0.197937
0.130951 5.35736e-002
1.17645e-002 0.363892
0.58551 0.26387
例子2代码：
f(x1,x2,x3,y1,y2,y3)= //函数定义
{
y1=x1*x1+x2*x2+x3*x3-1.0,
y2=2.0*x1*x1+x2*x2-4.0*x3,
y3=3.0*x1*x1-4.0*x2+x3*x3
};
!using["math","sys"];
mvar:
oo{
x=array(3),
x.SA, //设置初值为1,1,1
i=netn, //拟牛顿法解方程
x.outm(), //输出结果
i //返回迭代次数
};
结果：
0.785197 0.496611 0.369923
效率测试：
simwe的网友lin2009 的matlab代码：
clear all
clc
tic
k = zeros(5,5); % //生成5×5全0矩阵
% 循环计算以下程序段100000次：
for m = 1:100000
a = rand(5,7);
b = rand(7,5);%//生成5×7矩阵a,7×5矩阵b,用0～1之间的随机数初始化
k = k + a * b + a(1:5, 2:6) * b(2:6, 1:5) - a(:, 7) * b(3, :);
end
k
toc
Forcal代码：
运行稍快的代码，比matlab约快10％吧？
!using["math","sys"];
mvar:
t0=clock(),
oo{k=zeros}, //生成5×5矩阵k，初始化为0
i=0,(i<1000 00).while{ //循环计算1000 00次
oo{
a=rand, b=rand, //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
k.oset //计算k=k+a*b+a.subg(0,4:1,5)*b.subg(1,5:0,4)-a.subg(neg:6)*b.subg(3:neg)
},
i++
},
k.outm(), //输出矩阵k，然后销毁k
/1000; //得到计算时间，秒
在我的电脑上运行时间为3.344秒。
比较好看些的代码，似乎也比matlab稍快吧？
!using["math","sys"];
(:t0,k,i,a,b)=
{
t0=clock(),
k=zeros,
i=0,(i<1000 00).while{
oo{
a=rand, b=rand,
k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)-a(neg:6)*b(3:neg)
},
i++
},
k.outm().delete(),
/1000
};
在我的电脑上运行时间为3.579秒。
该例子的理论结果是每个元素均为275000。
我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。
继续例子，大家看有什么问题吗？
!using["math"];
mvar:
oo{
ndgrid,
a1.outm,
a2.outm,
a3.outm,
a4.outm,
a=a1+a2+a3+a4,
a.outm,
Sum.outm.Sum[].outm.Sum[].outm.Sum[]
};
说明：
linspace(8,12,5)：生成一维数组，共5个元素8～12
a1.outm：输出**数组a1，连下标一起输出
Sum：设有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.
在matlab中，纯for循环速度最慢。而一半for循环+一半向量化的速度最快，Forcal中也是如此：
!using["math","sys"];
mvar:
(:p1,p2,p3,a,b)=
{
oo{
a=array.rand(),
b=array.rand(),
p1=array,
p2=array,
p3=array,
t0=clock(),
ndgrid(a,b,&A,&B),
p1.=A+B
},
printff{"\r\nndgrid: {1,r}",/1000},
lena=FCDLen(a),
lenb=FCDLen(b),
t0=clock(),
m = lenb-1, (m>=0).while{
oo{p2(m,neg) = a+rn},
m--
},
printff{"\r\nfor1: {1,r}",/1000},
t0=clock(),
m = lenb-1, (m>=0).while{
n = lena-1, (n>=0).while{
//p3(m,n) = a(n)+b(m), //用这句还要慢一些
A(p3,m,n) = A(a,n)+A(b,m),
n--
},
m--
},
printff{"\r\nfor2: {1,r}",/1000}
};
结果：
ndgrid: 3.2001e-002
for1: 1.4999e-002
for2: 1.86
一段程序的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"];
mvar:
t=clock(),
oo{
ndgrid,
Sum
};
/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)))));
sum["f",0,1,0.0011,1,2,0.0009];
结果：
1008606.64947441
0.719 //时间
3、while循环
mvar:
t=sys::clock();
s=0,x=0,
while{x<=1, //while循环算法；
y=1,
while{y<=2,
s=s+cos(1-sin(1.2*x^(y/2)+cos(1-sin(1.2*y^(x/2))))),
y=y+0.0009
},
x=x+0.0011
},
s;
/1000;
结果：
1008606.64947441
0.734 //时间
好深奥！~~~~ 本帖最后由 forcal 于 2010-10-8 21:10 编辑
好深奥！~~~~
qbist 发表于 2010-10-7 14:56 http://www.madio.net/static/image/common/back.gif
先了解一下，以备不时之需，有问题可以交流哦，呵呵。 回复 forcal 的帖子
嗯！！！ 改进了FcMath中的矩阵乘算法，不知与matlab还有多大差距，朋友们可帮助测一下。
以下是FcMath中的矩阵乘与徐士良算法库XSLSF（普通的C/C++算法）中的矩阵乘的效率比较：
1、FcMath中的矩阵乘
!using["math","sys"];
(:a,b,k,t0)=
oo{
a=rand, b=rand,
t0=clock(),
k=a*b, //矩阵乘
k.outm()
},
/1000;
结果：
238.447 247.837 247.065 248.105 247.058
244.123 249.925 247.553 243.981 250.016
236.387 252.025 245.651 248.866 248.866
2.219 秒
2、XSLSF（普通的C/C++算法）中的矩阵乘
!using["math","sys","XSLSF"];
(:a,b,k,t0)=
oo{
a=rand, b=rand, k=array,
t0=clock(),
rmul, //矩阵乘
k.outm()
},
/1000;
结果：
262.121 247.583 260.529 259.548 258.328
255.413 246.563 254.356 250.548 251.509
256.152 247.725 259.444 250.827 249.816
10.563 秒
页:
[1]