|
楼主 |
发表于 2005-1-19 22:42:48
|
显示全部楼层
//复系数全选主元高斯-约当消去法
//ar-n*n 存放方程组的系数矩阵的实部,返回时将被破坏
//ai-n*n 存放方程组的系数矩阵的虚部,返回时将被破坏
//br-n*m常数向量的实部,每列为一组,返回时存放m组解向量的实部
//bi-n*m常数向量的虚部,每列为一组,返回时存放m组解向量的虚部
//n-方程组的阶数
//m-方程组右端常数向量的个数
//返回0表示原方程组的系数矩阵奇异
int cdcgj(double ar[],double ai[],double br[],double bi[],int n,int m)
{
int *js,l,k,i,j,is,u,v;
double p,q,s,d;
js=malloc(n*sizeof(int));
for (k=0;k<=n-1;k++)
{
d=0.0;
for (i=k;i<=n-1;i++)
{
for (j=k;j<=n-1;j++)
{
u=i*n+j;
p=ar*ar+ai*ai;
if (p>d)
{
d=p;js[k]=j;is=i;
}
}
}
if (d+1.0==1.0)
{
free(js);
printf("err**fail\n");
return(0);
}
if (is!=k)
{
for (j=k;j<=n-1;j++)
{
u=k*n+j; v=is*n+j;
p=ar; ar=ar[v]; ar[v]=p;
p=ai; ai=ai[v]; ai[v]=p;
}
for (j=0;j<=m-1;j++)
{
u=k*m+j; v=is*m+j;
p=br; br=br[v]; br[v]=p;
p=bi; bi=bi[v]; bi[v]=p;
}
}
if (js[k]!=k)
{
for (i=0;i<=n-1;i++)
{
u=i*n+k; v=i*n+js[k];
p=ar; ar=ar[v]; ar[v]=p;
p=ai; ai=ai[v]; ai[v]=p;
}
}
v=k*n+k;
for (j=k+1;j<=n-1;j++)
{
u=k*n+j;
p=ar*ar[v]; q=-ai*ai[v];
s=(ar[v]-ai[v])*(ar+ai);
ar=(p-q)/d; ai=(s-p-q)/d;
}
for (j=0;j<=m-1;j++)
{
u=k*m+j;
p=br*ar[v]; q=-bi*ai[v];
s=(ar[v]-ai[v])*(br+bi);
br=(p-q)/d; bi=(s-p-q)/d;
}
for (i=0;i<=n-1;i++)
{
if (i!=k)
{
u=i*n+k;
for (j=k+1;j<=n-1;j++)
{
v=k*n+j; l=i*n+j;
p=ar*ar[v]; q=ai*ai[v];
s=(ar+ai)*(ar[v]+ai[v]);
ar[l]=ar[l]-p+q;
ai[l]=ai[l]-s+p+q;
}
for (j=0;j<=m-1;j++)
{
l=i*m+j; v=k*m+j;
p=ar*br[v]; q=ai*bi[v];
s=(ar+ai)*(br[v]+bi[v]);
br[l]=br[l]-p+q; bi[l]=bi[l]-s+p+q;
}
}
}
}
for (k=n-1;k>=0;k--)
{
if (js[k]!=k)
{
for (j=0;j<=m-1;j++)
{
u=k*m+j;v=js[k]*m+j;
p=br; br=br[v]; br[v]=p;
p=bi; bi=bi[v]; bi[v]=p;
}
}
}
free(js);
return(1);
}
|
|