这个高斯列主元消去法,哪里错了?
#include<stdio.h>#include<stdlib.h>
#include<math.h>
#define N 4
#define Error 0.00001
void Gauss(float p[N][N+1],float *p1);
void Solve(float p[N][N+1],float *p1,float *p2,float *p3,int rank,int a_rank);
void Print(float p[N][N+1]);
void Solve_2(float p[N][N+1],float *p1,int rank,float *p3);
int main(void)
{
int i,j;
float p[N][N+1],*p1=NULL; //用于调用高斯列主元消去法后，返回并打印方程唯一解
p1=(float*)calloc(N,sizeof(float));
if(p1==NULL)
{
printf("No enough memory(main_p)!\n");
exit(1);
}
for(i=0;i<N;i++) //输入增广矩阵
{
printf("Row %d,enter %d elements:\n",i+1,N+1);
for(j=0;j<N+1;j++)
{
scanf("%f",&p[i][j]);
}
}
Gauss(p,p1); //调用高斯列主元消去法
printf("\n");
free(p1);
system("pause");
return(0);
}
void Gauss(float p[N][N+1],float *p1) //p1用于保存解
{
int i,j,k,index,revise=0,s,t,rank=0,a_rank=0; //index是某一列主元的下标，revise起到修正主元下标的作用
float *p2=NULL,*p3=NULL,Coeff; //p2用于换行以及打印基础解系，p3用于保存主元及其所在的列
p2=(float*)calloc(N,sizeof(float));
p3=(float*)calloc(N,sizeof(float));
if(p2==NULL)
{
printf("No enough memory(Gauss_p2)!\n");
exit(1);
}
if(p3==NULL)
{
printf("No enough memory(Gauss_p3)!\n");
exit(1);
}
for(j=0;j<N;j++) //表示第j列选主元
{
index=j-revise;
for(i=j+1-revise;i<N;i++) //找到第j列主元的下标位置
{
if(fabs(p[index][j])+Error<fabs(p[i][j]))
index=i;
}
for(k=j;k<N+1;k++) //第j-revise行与主元所在的行交换
{
p2[k]=p[j-revise][k];
p[j-revise][k]=p[index][k];
p[index][k]=p2[k];
}
if(fabs(p[j-revise][j])<Error) //p[j-revise][j]是主元
{
revise=revise+1;
continue;
}
else
{
rank=rank+1; //计算系数矩阵的秩
p3[j]=p[j-revise][j]; //把主元保存在p3数组中
}
for(s=j+1-revise;s<N;s++) //消元
{
Coeff=p[s][j]/p[j-revise][j];
for(t=j;t<N+1;t++)
p[s][t]=p[s][t]-Coeff*p[j-revise][t];
}
}
/*printf("\nPrimary:\n"); //打印p3数组，显示主元和主元所在的列
for(i=0;i<N;i++)
printf("%f ",p3[i]);
printf("\n");*/
for(i=N;i>=0;i--) //计算增广矩阵的秩
{
if(fabs(p[i][N])>Error)
{
a_rank=i+1;
break;
}
}
Solve(p,p1,p2,p3,rank,a_rank);
free(p2);
free(p3);
}
void Solve(float p[N][N+1],float *p1,float *p2,float *p3,int rank,int a_rank)
{
int i,j,k;
float term,sum;
if(rank==N) //系数矩阵满秩 ，方程有唯一解
{
printf("\nAfter Gauss:\n"); //打印经过高斯列主元消去法后的增广矩阵
Print(p);
for(i=N-1;i>=0;i--)
{
sum=0;
for(j=N-1;j>i;j--)
{
term=p[i][j]*p1[j];
sum=sum+term;
}
p1[i]=(p[i][N]-sum)/p3[i];
}
printf("\nSolution:"); //打印方程的唯一解
for(i=0;i<N;i++)
printf("%f ",p1[i]);
printf("\n");
}
else if(rank>=a_rank && rank<N) //系数矩阵不满秩，系数矩阵的秩等于增广矩阵的秩，方程有无穷多解
{
printf("\nAfter Gauss:\n"); //打印经过高斯列主元消去法后的增广矩阵
Print(p);
Solve_2(p,p1,rank-1,p3);
for(i=N-1;i>=0;i--)
{
sum=0;
for(j=N-1;j>i;j--)
{
term=p[i][j]*p1[j];
sum=sum+term;
}
if(fabs(p3[i])>Error)
p1[i]=(p[i][N]-sum)/p3[i];
}
printf("\nParticular solution:\n");
for(i=0;i<N;i++)
printf("%f ",p1[i]);
for(i=0;i<N;i++)
p[i][N]=0;
for(k=N-1;k>=0;k--)
{
if(fabs(p3[k])<Error)
{
p1[k]=1;
for(i=N-1;i>=0;i--)
{
sum=0;
for(j=N-1;j>i;j--)
{
term=p[i][j]*p1[j];
sum=sum+term;
}
if(fabs(p3[i])>Error)
p1[i]=(p[i][N]-sum)/p3[i];
}
printf("\nGeneral solution:\n");
for(i=0;i<N;i++)
printf("%f ",p1[i]);
p1[k]=0;
}
}
}
else if(rank<a_rank) //系数矩阵的秩不等于增广矩阵的秩，方程无解
{
printf("\nAfter Gauss:\n"); //打印经过高斯列主元消去法后的增广矩阵
Print(p);
printf("\nNo solution!\n");
}
}
void Print(float p[N][N+1])
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<N+1;j++)
printf("%f ",p[i][j]);
printf("\n");
}
}
void Solve_2(float p[N][N+1],float *p1,int rank,float *p3) //p1和p3就是Solve中的p1和p3
{
int i,j,k,s,t;
float *p2=NULL,*p4=NULL;
p2=(float*)calloc(N,sizeof(float)); //用来换行
p4=(float*)calloc(N,sizeof(float)); //用来保存p3
if(p2==NULL)
{
printf("No enough memory(Solve_p2)!\n");
exit(1);
}
if(p4==NULL)
{
printf("No enough memory(Solve_p4)!\n");
exit(1);
}
for(i=0;i<N;i++)
p4[i]=p3[i];
for(i=rank;i>=0;i--)
{
for(j=N-1;j>=0;j--)
{
if(fabs(p4[j])>Error)
{
for(k=j;k<N+1;k++)
{
p2[k]=p[j][k];
p[j][k]=p[i][k];
p[i][k]=p2[k];
}
p4[j]=0;
break;
}
}
}
free(p2);
free(p4);
}
这个是用Dev C++写的，同样一组数据，有时候可以运行，有时候不可以。。。
