极限测试之Matlab与Forcal真实演练
首先要说明的是本测试系列，除了比较编译效率外，其余所有运行效率的比较都是将matlab的首次运行排出在外，因matlab程序首次运行效率较低。理论上，Forcal程序任意次运行效率都是一样的。不过话又说回来，任意程序包含的函数，有些是需要多次运行的，而有些仅运行一次，甚至一次都不运行，故matlab函数首次运行效率较低应该是一个缺点。但如果说，matlab函数首次运行会对函数进行优化，以后运行效率会显著提高，则matlab函数首次运行效率较低就成了一个优点。＝＝＝＝＝＝＝＝＝＝＝＝＝
本次演练要用matlab和Forcal实现两个实用函数并进行测试。当然，对脚本来说，这些函数用C/C++或Fortran来实现应是脚本的最佳选择。
＝＝＝＝＝＝＝＝＝＝＝＝＝
1、求解实系数方程组的全选主元高斯消去法：包含大量数组元素存取操作
C/C++代码：#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include "time.h"
#include "math.h"
int agaus(double *a,double *b,int n)
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=new int;
l=1;
for (k=0;k<=n-2;k++)
{
d=0.0;
for (i=k;i<=n-1;i++)
{
for (j=k;j<=n-1;j++)
{
t=fabs(a);
if (t>d) { d=t; js=j; is=i;}
}
}
if (d+1.0==1.0)
{
l=0;
}
else
{
if (js!=k)
{
for (i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js;
t=a; a=a; a=t;
}
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j;
t=a; a=a; a=t;
}
t=b; b=b; b=t;
}
}
if (l==0)
{
delete[] js; printf("fail\n");
return(0);
}
d=a;
for (j=k+1;j<=n-1;j++)
{
p=k*n+j; a=a/d;
}
b=b/d;
for (i=k+1;i<=n-1;i++)
{
for (j=k+1;j<=n-1;j++)
{
p=i*n+j;
a=a-a*a;
}
b=b-a*b;
}
}
d=a[(n-1)*n+n-1];
if (fabs(d)+1.0==1.0)
{
delete[] js; printf("fail\n");
return(0);
}
b=b/d;
for (i=n-2;i>=0;i--)
{
t=0.0;
for (j=i+1;j<=n-1;j++)
{
t=t+a*b;
}
b=b-t;
}
js=n-1;
for (k=n-1;k>=0;k--)
{
if (js!=k)
{
t=b; b=b]; b]=t;
}
}
delete[] js;
return(1);
}
int main(int argc, char *argv[])
{
int i,j,k;
double a=
{ {0.2368,0.2471,0.2568,1.2671},
{0.1968,0.2071,1.2168,0.2271},
{0.1581,1.1675,0.1768,0.1871},
{1.1161,0.1254,0.1397,0.1490} };
double b={1.8471,1.7471,1.6471,1.5471};
double aa,bb;
clock_t tm;
tm=clock();
for(i=0;i<10000;i++)
{
for(j=0;j<4;j++)
{
for(k=0;k<4;k++)
{
aa=a;
}
}
for(j=0;j<4;j++)
{
bb=b;
}
agaus((double *)aa,bb,4);
}
printf("循环 %d 次, 耗时 %d 毫秒。\n", i,(clock()-tm));
for (i=0;i<=3;i++)
{
printf("x(%d)=%e\n",i,bb);
}
}结果：
循环 10000 次, 耗时 31 毫秒。
x(0)=1.040577e+000
x(1)=9.870508e-001
x(2)=9.350403e-001
x(3)=8.812823e-001
－－－－－－－－－
matlab 2009a代码：%file agaus.m
function c=agaus(a,b,n)
js=linspace(0,0,n);
l=1;
for k=1:n-1
d=0.0;
for i=k:n
for j=k:n
t=abs(a(i,j));
if (t>d)
d=t; js(k)=j; is=i;
end
end
end
if d+1.0==1.0
l=0;
else
if js(k)~=k
for i=1:n
t=a(i,k); a(i,k)=a(i,js(k)); a(i,js(k))=t;
end
end
if is~=k
for j=k:n
t=a(k,j); a(k,j)=a(is,j); a(is,j)=t;
end
t=b(k); b(k)=b(is); b(is)=t;
end
end
if l==0
printf('fail\n');
c=[];
return;
end
d=a(k,k);
for j=k+1:n
a(k,j)=a(k,j)/d;
end
b(k)=b(k)/d;
for i=k+1:n
for j=k+1:n
a(i,j)=a(i,j)-a(i,k)*a(k,j);
end
b(i)=b(i)-a(i,k)*b(k);
end
end
d=a(n,n);
if abs(d)+1.0==1.0
printf('fail\n');
c=[];
return;
end
b(n)=b(n)/d;
for i=n-1:-1:1
t=0.0;
for j=i+1:n
t=t+a(i,j)*b(j);
end
b(i)=b(i)-t;
end
js(n)=n;
for k=n:-1:1
if js(k)~=k
t=b(k); b(k)=b(js(k)); b(js(k))=t;
end
end
c=b;
return;
end
a=[0.2368,0.2471,0.2568,1.2671;
0.1968,0.2071,1.2168,0.2271;
0.1581,1.1675,0.1768,0.1871;
1.1161,0.1254,0.1397,0.1490] ;
b=[ 1.8471,1.7471,1.6471,1.5471];
tic
for i=1:10000
c=agaus(a,b,4);
end
c
toc
c =
1.0406 0.9871 0.9350 0.8813
Elapsed time is 0.762713 seconds.－－－－－－－－－－
Forcal代码：!using["math","sys"];
agaus(a,b,n : js,l,k,i,j,is, d,t)=
{
oo{ js=array(n)},
l=1, k=0,
while{ k<n-1,
d=0.0, i=k,
while{ i<n,
j=k, while{j<n,
t=abs(a),
if{t>d, d=t, js=j, is=i},
j++
},
i++
},
which{ d+1.0==1.0, l=0,
{ if{ (js!=k),
i=0, while{i<n,
t=a, a=a], a]=t,
i++
}
},
if{ (is!=k),
j=k, while{j<n,
t=a, a=a, a=t,
j++
},
t=b, b=b, b=t
}
}
},
if{ (l==0),
printff("fail\r\n"),
return(0)
},
d=a,
j=k+1, while {j<n, a=a/d, j++},
b=b/d,
i=k+1, while {i<n,
j=k+1, while{j<n,
a=a-a*a,
j++
},
b=b-a*b,
i++
},
k++
},
d=a[(n-1),n-1],
if{ abs(d)+1.0==1.0,
printff("fail\r\n"),
return(0)
},
b=b/d,
i=n-2, while{i>=0,
t=0.0,
j=i+1, while{j<n, t=t+a*b, j++},
b=b-t,
i--
},
js=n-1,
k=n-1, while{k>=0,
if{(js!=k), t=b, b=b], b]=t},
k--
},
return(1)
};
main(:i,a,b,aa,bb,t0)=
{
oo{a=arrayinit{2,4,4 :
0.2368,0.2471,0.2568,1.2671,
0.1968,0.2071,1.2168,0.2271,
0.1581,1.1675,0.1768,0.1871,
1.1161,0.1254,0.1397,0.1490},
b=arrayinit{1,4 : 1.8471,1.7471,1.6471,1.5471},
aa=array, bb=array
},
t0=clock(),
i=0, while{i<10000, aa.=a, bb.=b, agaus(aa,bb,4), i++},
outm,
/1000
};结果：
1.04058 0.987051 0.93504 0.881282
2.125
Forcal用函数sys::A()对数组元素进行存取：!using["math","sys"];
agaus(a,b,n : js,l,k,i,j,is, d,t)=
{
oo{ js=array(n)},
l=1, k=0,
while{ k<n-1,
d=0.0, i=k,
while{ i<n,
j=k, while{j<n,
t=abs(A),
if{t>d, d=t, A=j, is=i},
j++
},
i++
},
which{ d+1.0==1.0, l=0,
{ if{ (A!=k),
i=0, while{i<n,
t=A, A=A], A]=t,
i++
}
},
if{ (is!=k),
j=k, while{j<n,
t=A, A=A, A=t,
j++
},
t=A, A=A, A=t
}
}
},
if{ (l==0),
printff("fail\r\n"),
return(0)
},
d=A,
j=k+1, while {j<n, A=A/d, j++},
A=A/d,
i=k+1, while {i<n,
j=k+1, while{j<n,
A=A-A*A,
j++
},
A=A-A*A,
i++
},
k++
},
d=A,
if{ abs(d)+1.0==1.0,
printff("fail\r\n"),
return(0)
},
A=A/d,
i=n-2, while{i>=0,
t=0.0,
j=i+1, while{j<n, t=t+A*A, j++},
A=A-t,
i--
},
A=n-1,
k=n-1, while{k>=0,
if{(A!=k), t=A, A=A], A]=t},
k--
},
return(1)
};
main(:i,a,b,aa,bb,t0)=
{
oo{a=arrayinit{2,4,4 :
0.2368,0.2471,0.2568,1.2671,
0.1968,0.2071,1.2168,0.2271,
0.1581,1.1675,0.1768,0.1871,
1.1161,0.1254,0.1397,0.1490},
b=arrayinit{1,4 : 1.8471,1.7471,1.6471,1.5471},
aa=array, bb=array
},
t0=clock(),
i=0, while{i<10000, aa.=a, bb.=b, agaus(aa,bb,4), i++},
outm,
/1000
};结果：
1.04058 0.987051 0.93504 0.881282
1.454
－－－－－－－－－－
可以看出C/C++、matlab、Forcal耗时之比为 1 ：25：68 （Forcal不使用函数sys::A()）。
可以看出C/C++、matlab、Forcal耗时之比为 1 ：25：47 （Forcal使用函数sys::A()）。
本例Forcal耗时较长的原因在于本例程序含有大量的数组元素存取操作。 2、变步长辛卜生二重求积法：没有数组元素操作
C/C++代码：#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include "time.h"
#include "math.h"
double simp1(double x,double eps);
void fsim2s(double x,double y[]);
double fsim2f(double x,double y);
double fsim2(double a,double b,double eps)
{
int n,j;
double h,d,s1,s2,t1,x,t2,g,s,s0,ep;
n=1; h=0.5*(b-a);
d=fabs((b-a)*1.0e-06);
s1=simp1(a,eps); s2=simp1(b,eps);
t1=h*(s1+s2);
s0=1.0e+35; ep=1.0+eps;
while (((ep>=eps)&&(fabs(h)>d))||(n<16))
{
x=a-h; t2=0.5*t1;
for (j=1;j<=n;j++)
{
x=x+2.0*h;
g=simp1(x,eps);
t2=t2+h*g;
}
s=(4.0*t2-t1)/3.0;
ep=fabs(s-s0)/(1.0+fabs(s));
n=n+n; s0=s; t1=t2; h=h*0.5;
}
return(s);
}
double simp1(double x,double eps)
{
int n,i;
double y,h,d,t1,yy,t2,g,ep,g0;
n=1;
fsim2s(x,y);
h=0.5*(y-y);
d=fabs(h*2.0e-06);
t1=h*(fsim2f(x,y)+fsim2f(x,y));
ep=1.0+eps; g0=1.0e+35;
while (((ep>=eps)&&(fabs(h)>d))||(n<16))
{
yy=y-h;
t2=0.5*t1;
for (i=1;i<=n;i++)
{
yy=yy+2.0*h;
t2=t2+h*fsim2f(x,yy);
}
g=(4.0*t2-t1)/3.0;
ep=fabs(g-g0)/(1.0+fabs(g));
n=n+n; g0=g; t1=t2; h=0.5*h;
}
return(g);
}
void fsim2s(double x,double y[])
{
y=-sqrt(1.0-x*x);
y=-y;
}
double fsim2f(double x,double y)
{
return exp(x*x+y*y);
}
int main(int argc, char *argv[])
{
int i;
double a,b,eps,s;
clock_t tm;
a=0.0; b=1.0; eps=0.0001;
tm=clock();
for(i=0;i<100;i++)
{
s=fsim2(a,b,eps);
}
printf("s=%e , 耗时 %d 毫秒。\n", s, (clock()-tm));
}结果：
s=2.698925e+000 , 耗时 78 毫秒。
－－－－－－－
matlab代码：%file fsim2.m
function s=fsim2(a,b,eps)
n=1; h=0.5*(b-a);
d=abs((b-a)*1.0e-06);
s1=simp1(a,eps); s2=simp1(b,eps);
t1=h*(s1+s2);
s0=1.0e+35; ep=1.0+eps;
while ((ep>=eps)&&(abs(h)>d))||(n<16),
x=a-h; t2=0.5*t1;
for j=1:n
x=x+2.0*h;
g=simp1(x,eps);
t2=t2+h*g;
end
s=(4.0*t2-t1)/3.0;
ep=abs(s-s0)/(1.0+abs(s));
n=n+n; s0=s; t1=t2; h=h*0.5;
end
end
function g=simp1(x,eps)
n=1;
=f2s(x);
h=0.5*(y1-y0);
d=abs(h*2.0e-06);
t1=h*(f2f(x,y0)+f2f(x,y1));
ep=1.0+eps; g0=1.0e+35;
while (((ep>=eps)&&(abs(h)>d))||(n<16))
yy=y0-h;
t2=0.5*t1;
for i=1:n
yy=yy+2.0*h;
t2=t2+h*f2f(x,yy);
end
g=(4.0*t2-t1)/3.0;
ep=abs(g-g0)/(1.0+abs(g));
n=n+n; g0=g; t1=t2; h=0.5*h;
end
end
%file f2s.m
function =f2s(x)
y0=-sqrt(1.0-x*x);
y1=-y0;
end
%file f2f.m
function c=f2f(x,y)
c=exp(x*x+y*y);
end
%%%%%%%%%%%%%
>> tic
for i=1:100
a=fsim2(0,1,0.0001);
end
a
toc
a =
2.6989
Elapsed time is 0.995575 seconds.－－－－－－－
Forcal代码：fsim2s(x,y0,y1)=
{
y0=-sqrt(1.0-x*x),
y1=-y0
};
fsim2f(x,y)=exp(x*x+y*y);
//////////////////
simp1(x,eps : n,i,y0,y1,h,d,t1,yy,t2,g,ep,g0)=
{
n=1,
fsim2s(x,&y0,&y1),
h=0.5*(y1-y0),
d=abs(h*2.0e-06),
t1=h*(fsim2f(x,y0)+fsim2f(x,y1)),
ep=1.0+eps, g0=1.0e+35,
while {((ep>=eps)&(abs(h)>d))|(n<16),
yy=y0-h,
t2=0.5*t1,
i=1, while{i<=n,
yy=yy+2.0*h,
t2=t2+h*fsim2f(x,yy),
i++
},
g=(4.0*t2-t1)/3.0,
ep=abs(g-g0)/(1.0+abs(g)),
n=n+n, g0=g, t1=t2, h=0.5*h
},
g
};
fsim2(a,b,eps : n,j,h,d,s1,s2,t1,x,t2,g,s,s0,ep)=
{
n=1, h=0.5*(b-a),
d=abs((b-a)*1.0e-06),
s1=simp1(a,eps), s2=simp1(b,eps),
t1=h*(s1+s2),
s0=1.0e+35, ep=1.0+eps,
while {((ep>=eps)&(abs(h)>d))|(n<16),
x=a-h, t2=0.5*t1,
j=1, while{j<=n,
x=x+2.0*h,
g=simp1(x,eps),
t2=t2+h*g,
j++
},
s=(4.0*t2-t1)/3.0,
ep=abs(s-s0)/(1.0+abs(s)),
n=n+n, s0=s, t1=t2, h=h*0.5
},
s
};
//////////////////
mvar:
t0=sys::clock(),
i=0, while{i<100, a=fsim2(0,1,0.0001), i++}, a;
/1000;结果：
2.698925000624303
0.328
－－－－－－－－－
本例C/C++、matlab、Forcal运行耗时之比为 1：12.7：4.2 。
本例matlab慢的原因应该在于有大量的函数调用。另外，matlab要求每一个函数必须存为磁盘文件真的很麻烦。
本例还说明，对C/C++和Forcal脚本来说，C/C++的一条指令，Forcal平均大约用4.2条指令进行解释。 3、变步长辛卜生二重求积法，写成通用的函数：没有数组元素操作
注意函数fsim2中增加了两个参数，函数句柄fsim2s用于计算二重积分时内层的上下限，函数句柄fsim2f用于计算积分函数的值。
不再给出C/C++代码，因其效率不会发生变化。
Matlab代码：%file fsim2.m
function s=fsim2(a,b,eps,fsim2s,fsim2f)
n=1; h=0.5*(b-a);
d=abs((b-a)*1.0e-06);
s1=simp1(a,eps,fsim2s,fsim2f); s2=simp1(b,eps,fsim2s,fsim2f);
t1=h*(s1+s2);
s0=1.0e+35; ep=1.0+eps;
while ((ep>=eps)&&(abs(h)>d))||(n<16),
x=a-h; t2=0.5*t1;
for j=1:n
x=x+2.0*h;
g=simp1(x,eps,fsim2s,fsim2f);
t2=t2+h*g;
end
s=(4.0*t2-t1)/3.0;
ep=abs(s-s0)/(1.0+abs(s));
n=n+n; s0=s; t1=t2; h=h*0.5;
end
end
function g=simp1(x,eps,fsim2s,fsim2f)
n=1;
=fsim2s(x);
h=0.5*(y1-y0);
d=abs(h*2.0e-06);
t1=h*(fsim2f(x,y0)+fsim2f(x,y1));
ep=1.0+eps; g0=1.0e+35;
while (((ep>=eps)&&(abs(h)>d))||(n<16))
yy=y0-h;
t2=0.5*t1;
for i=1:n
yy=yy+2.0*h;
t2=t2+h*fsim2f(x,yy);
end
g=(4.0*t2-t1)/3.0;
ep=abs(g-g0)/(1.0+abs(g));
n=n+n; g0=g; t1=t2; h=0.5*h;
end
end
%file f2s.m
function =f2s(x)
y0=-sqrt(1.0-x*x);
y1=-y0;
end
%file f2f.m
function c=f2f(x,y)
c=exp(x*x+y*y);
end
%%%%%%%%%%%%%%%%
>> tic
for i=1:100
a=fsim2(0,1,0.0001,@f2s,@f2f);
end
a
toc
a =
2.6989
Elapsed time is 1.267014 seconds.－－－－－－－－
Forcal代码：simp1(x,eps,fsim2s,fsim2f : n,i,y0,y1,h,d,t1,yy,t2,g,ep,g0)=
{
n=1,
fsim2s(x,&y0,&y1),
h=0.5*(y1-y0),
d=abs(h*2.0e-06),
t1=h*(fsim2f(x,y0)+fsim2f(x,y1)),
ep=1.0+eps, g0=1.0e+35,
while {((ep>=eps)&(abs(h)>d))|(n<16),
yy=y0-h,
t2=0.5*t1,
i=1, while{i<=n,
yy=yy+2.0*h,
t2=t2+h*fsim2f(x,yy),
i++
},
g=(4.0*t2-t1)/3.0,
ep=abs(g-g0)/(1.0+abs(g)),
n=n+n, g0=g, t1=t2, h=0.5*h
},
g
};
fsim2(a,b,eps,fsim2s,fsim2f : n,j,h,d,s1,s2,t1,x,t2,g,s,s0,ep)=
{
n=1, h=0.5*(b-a),
d=abs((b-a)*1.0e-06),
s1=simp1(a,eps,fsim2s,fsim2f), s2=simp1(b,eps,fsim2s,fsim2f),
t1=h*(s1+s2),
s0=1.0e+35, ep=1.0+eps,
while {((ep>=eps)&(abs(h)>d))|(n<16),
x=a-h, t2=0.5*t1,
j=1, while{j<=n,
x=x+2.0*h,
g=simp1(x,eps,fsim2s,fsim2f),
t2=t2+h*g,
j++
},
s=(4.0*t2-t1)/3.0,
ep=abs(s-s0)/(1.0+abs(s)),
n=n+n, s0=s, t1=t2, h=h*0.5
},
s
};
//////////////////
f2s(x,y0,y1)=
{
y0=-sqrt(1.0-x*x),
y1=-y0
};
f2f(x,y)=exp(x*x+y*y);
mvar:
t0=sys::clock(),
i=0, while{i<100, a=fsim2(0,1,0.0001,HFor("f2s"),HFor("f2f")), i++}, a;
/1000;结果：
2.698925000624303
0.844
－－－－－－－－
本例matlab与Forcal耗时之比为1.267014 ：0.844。Forcal仍有优势。
本例Forcal耗时增加的原因：在函数fsim2及simp1中要动态查找函数句柄fsim2s,fsim2f，并验证其是否有效，故效率下降了。 百年不遇的好帖子，不得不顶 确实很不错 好深奥哦，肿么办？
页:
[1]