gauss.cpp

实现功能：顺序高斯消去法、列主元素消去法、全主元素消去法

//该程序由goldsen编写，仅供学习交流之用[e-mail]
//参考课本第二章例一，实现功能：顺序高斯消去法、列主元素消去法、全主元素消去法
//水平有限，不妥之处请批评指正，转载请注明出处:)
//////////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
////定义外部变量///////////////////
	double A[10][10];
	double B[10];
	double X[10];
    double b1=0.6781,b2=12.1,b3=981;
	signed int n=3;//选择矩阵的阶次
	
///////////////////////////////////
void input_matrix()
{

    
	signed int c;//选择矩阵输入方式
	static signed int i=1;//输出行
	signed int j=1;//输出列
	double a11=0.167,a12=0.01,a13=0.012,a21=5.91,a22=0.8334,a23=1,a31=4.2,a32=1200,a33=3200;//默认三阶矩阵
	
	double a=0;
	///////////////初始化矩阵空间,矩阵初始大小取决于外部变量定义,预设最大矩阵10,遇到不同情况可作更改
	for (i=1;i<=10;i++)
	{
		for (j=1;j<=10;j++)
		{
			A[i-1][j-1]=0;
		}
		B[i-1]=0;
		X[i-1]=0;
	}
	printf("Please insert your matrix.(0=default,1=insert)\n");
	scanf("%d",&c);
	if (0==c)
	{
		printf("you have chosen the default marix.\n");
		printf("The matrix is:\n");

		A[0][0]=a11;
		A[0][1]=a12;
		A[0][2]=a13;
		A[1][0]=a21;
		A[1][1]=a22;
		A[1][2]=a23;
		A[2][0]=a31;
		A[2][1]=a32;
		A[2][2]=a33;
		B[0]=b1;
		B[1]=b2;
		B[2]=b3;
		
		for (i=1;i<=3;i++)
		{
			for (j=1;j<=3;j++)
			{
				printf("%f \t",A[i-1][j-1]);
			}
			printf("\n");
		}
		printf("vector b is:\n%f\n%f\n%f\n\n",b1,b2,b3);
	} 
	else if (1==c)
	{
		printf("Please insert rank(阶次) n.\n");
		scanf("%d",&n);
		printf("Please insert %d elements in your matrix.\n",n*n);
		for (i=1;i<=n;i++)
		{
			for (j=1;j<=n;j++)
			{
				scanf("%lf",&a);
				A[i-1][j-1]=a;
			}
		}
		printf("The matrix is:\n");
		for (i=1;i<=n;i++)
		{
			for (j=1;j<=n;j++)
			{
				printf("%f \t",A[i-1][j-1]);
			}
			printf("\n");
		}
		printf("Please insert %d elements in vector b.\n",n);
		for (i=1;i<=n;i++)
		{
		    scanf("%lf",&a);
			B[i-1]=a;
			
		}
		printf("vector b is:\n");
		for (i=1;i<=n;i++)
		{
			printf("%f\n",B[i-1]);
		}
		printf("\n");
	}
	else
	{
		printf("you insert a wrong number.\n");
	}

}

double abs_double(double abs_in)/////////////double型取绝对值
{
	double abs_out;
	if (abs_in<=0)
	{
		abs_out=-abs_in;///////////////////////////////////取绝对值////////////
	} 
	else
	{
		abs_out=abs_in;
	}
	return abs_out;
	
}

void Gauss_shunxu()///////////////////////////////////////////////////////顺序高斯消元法
{
	signed int k;//消元次数
	signed int i;//行控制
	signed int j;//列控制
	double flag;
	double e=0.0001;//控制精度
	double m=0;//消元因子
	double Sum=0;//回代求和
	
/////////消元过程
	for (k=0;k<n-1;k++)
	{

		flag=abs_double(A[k][k]);/////////////////////////////////////////////////判断Akk的值，使之非零。
		
		if (flag<=e)
		{
			printf("A[%d][%d] is zero,wrong!",k+1,k+1);
			return;	
		}
		////////////////////////////////////////////////////////////
		for (i=k;i<=n-1;i++)
		{

			m=-A[i+1][k]/A[k][k];
			for (j=k;j<=n-1;j++)
			{
				A[i+1][j]=A[i+1][j]+m*A[k][j];
			}
			B[i+1]=B[i+1]+m*B[k];
		}
		printf("step%d:\n",k+1);
		for (i=1;i<=n;i++)
		{
			for (j=1;j<=n;j++)
			{
				printf("%f \t",A[i-1][j-1]);
			}
			printf("\n");
		}
		printf("vector b is:\n");
		for (i=1;i<=n;i++)
		{
			printf("%f\n",B[i-1]);
		}
		printf("\n");
	}

///////////////////////////////////////////////回代过程
	flag=abs_double(A[n-1][n-1]);

	if (flag<=e)//////////////////////////////////////////////////////////////////////////判断Ann的值，使之非零。
	{
		printf("A[%d][%d] is zero,wrong!\n",n-1,n-1);
		return;	
	}

	//////////////////////////////////////////////////////////////////////////
	X[n-1]=B[n-1]/A[n-1][n-1];
	for (k=n-1;k>=1;k--)
	{
		for (j=n-1;j>=k;j--)
		{
			Sum+=A[k-1][j]*X[j];
		}
		X[k-1]=(B[k-1]-Sum)/A[k-1][k-1];
		Sum=0;
	}
	printf("The %d roots is:\n",n);
	for (i=1;i<=n;i++)
	{
		printf("%f\n",X[i-1]);
	}


}

void Gauss_liezhu()//////////////////////////////////////////////////////列主元素消去法
{
	signed int k;//消元次数
	signed int i;//行控制
	signed int j;//列控制
	double flag;
	double e=0.00001;//控制精度
	double m=0;//消元因子
	double Sum=0;//回代求和
	double max_A;
	double t_A;
	signed int ik;
	
/////////消元过程
	for (k=0;k<n-1;k++)
	{
		//////////////////////////////////////////////////////////////////////////交换k行和ik行
		
		for (i=k,max_A=A[k][k],ik=k;i<=n-1;i++)
		{
			if (abs_double(max_A)<abs_double(A[i][k]))
			{
				max_A=A[i][k];
				ik=i;
			}

		}
		for (j=0;j<=n-1;j++)
		{
			t_A=A[k][j];
			A[k][j]=A[ik][j];
			A[ik][j]=t_A;
		}

		t_A=B[k];
		B[k]=B[ik];
		B[ik]=t_A;

		

        //////////////////////////////////////////////////////////////////////////

		flag=abs_double(A[k][k]);/////////////////////////////////////////////////判断Akk的值，使之非零。
		
		if (flag<=e)
		{
			printf("A[%d][%d] is zero,wrong!",k+1,k+1);
			return;	
		}
		////////////////////////////////////////////////////////////
		for (i=k;i<=n-1;i++)
		{

			m=-A[i+1][k]/A[k][k];
			for (j=k;j<=n-1;j++)
			{
				A[i+1][j]=A[i+1][j]+m*A[k][j];
			}
			B[i+1]=B[i+1]+m*B[k];
		}
		printf("step%d:\n",k+1);
		for (i=1;i<=n;i++)
		{
			for (j=1;j<=n;j++)
			{
				printf("%f \t",A[i-1][j-1]);
			}
			printf("\n");
		}
		printf("vector b is:\n");
		for (i=1;i<=n;i++)
		{
			printf("%f\n",B[i-1]);
		}
		printf("\n");
	}

///////////////////////////////////////////////回代过程
	flag=abs_double(A[n-1][n-1]);

	if (flag<=e)//////////////////////////////////////////////////////////////////////////判断Ann的值，使之非零。
	{
		printf("A[%d][%d] is zero,wrong!\n",n-1,n-1);
		return;	
	}

	//////////////////////////////////////////////////////////////////////////
	X[n-1]=B[n-1]/A[n-1][n-1];
	for (k=n-1;k>=1;k--)
	{
		for (j=n-1;j>=k;j--)
		{
			Sum+=A[k-1][j]*X[j];
		}
		X[k-1]=(B[k-1]-Sum)/A[k-1][k-1];
		Sum=0;
	}
	printf("The %d roots is:\n",n);
	for (i=1;i<=n;i++)
	{
		printf("%f\n",X[i-1]);
	}


}
void Gauss_quanzhu()//////////////////////////////////////////////////////全主元素消去法
{
	signed int k;//消元次数
	signed int i;//行控制
	signed int j;//列控制
	double flag;
	double e=0.00001;//控制精度
	double m=0;//消元因子
	double Sum=0;//回代求和
	double max_A;
	double t_A;
	signed int ik,jk;
	signed int order[10];//初始化矩阵空间,记录全主元的变量顺序,预设最大矩阵10,遇到不同情况可作更改
	for (i=0;i<=n-1;i++)
	{
		order[i]=i;
	}
	
/////////消元过程
	for (k=0;k<n-1;k++)
	{
		
//////////////////////////////////////////////////////////////////////////查找k行k列之后元素绝对值最大数
		for (i=k,max_A=A[k][k],jk=k,ik=k;i<=n-1;i++)
		{
			for (j=k;j<=n-1;j++)
			{
				if (abs_double(max_A)<abs_double(A[i][j]))
				{
					ik=i;
					jk=j;
				}
			}
		}
		for (j=0;j<=n-1;j++)////交换行
		{
			t_A=A[k][j];
			A[k][j]=A[ik][j];
			A[ik][j]=t_A;
		}
		
		t_A=B[k];////交换b
		B[k]=B[ik];
		B[ik]=t_A;
		
		for (i=0;i<=n-1;i++)//交换列
		{
			t_A=A[i][k];
			A[i][k]=A[i][jk];
			A[i][jk]=t_A;
		}

		i=order[jk];//记录列变序号
		order[jk]=order[k];
		order[k]=i;

        //////////////////////////////////////////////////////////////////////////

		flag=abs_double(A[k][k]);/////////////////////////////////////////////////判断Akk的值，使之非零。
		
		if (flag<=e)
		{
			printf("A[%d][%d] is zero,wrong!",k+1,k+1);
			return;	
		}
		////////////////////////////////////////////////////////////
		for (i=k;i<=n-1;i++)
		{

			m=-A[i+1][k]/A[k][k];
			for (j=k;j<=n-1;j++)
			{
				A[i+1][j]=A[i+1][j]+m*A[k][j];
			}
			B[i+1]=B[i+1]+m*B[k];
		}
		printf("step%d:\n",k+1);
		for (i=1;i<=n;i++)
		{
			for (j=1;j<=n;j++)
			{
				printf("%f \t",A[i-1][j-1]);
			}
			printf("\n");
		}
		printf("vector b is:\n");
		for (i=1;i<=n;i++)
		{
			printf("%f\n",B[i-1]);
		}
		printf("\n");
	}

///////////////////////////////////////////////回代过程
	flag=abs_double(A[n-1][n-1]);

	if (flag<=e)//////////////////////////////////////////////////////////////////////////判断Ann的值，使之非零。
	{
		printf("A[%d][%d] is zero,wrong!\n",n-1,n-1);
		return;	
	}

	//////////////////////////////////////////////////////////////////////////
	X[n-1]=B[n-1]/A[n-1][n-1];
	for (k=n-1;k>=1;k--)
	{
		for (j=n-1;j>=k;j--)
		{
			Sum+=A[k-1][j]*X[j];
		}
		X[k-1]=(B[k-1]-Sum)/A[k-1][k-1];
		Sum=0;
	}
	printf("The %d roots is:\n",n);
	for (i=1;i<=n;i++)
	{
		printf("X[%d]=%f\n",order[i-1]+1,X[i-1]);
	}


}
void main()
{
	signed int c;
	input_matrix();
	printf("plsease choose the method for guass elimination.(0=顺序Guass,1=列主元素Guass,2=全主元Guass)\n");
	scanf("%d",&c);
	if (0==c)
	{
		Gauss_shunxu();
		return;
	}
	else if (1==c)
	{
		Gauss_liezhu();
		return;
	}
	else if (2==c)
	{
		Gauss_quanzhu();
		return;
	}
}