列主元消元法

#include<iostream>
#include<cmath>
using namespace std;
#define N 4
void display(double(*a)[N], double *b)
{
	for(int i=0;i<N;i++)
	{
		for(int j=0;j<N;j++)
		{
			cout<<a[i][j]<<" ";
		}
		cout<<b[i]<<endl;
	}
}
int main()
{
	double a[N][N]={
		{1.1348, 3.8326, 1.1651, 3.4017},
		{0.5301, 1.7875, 2.5330, 1.5435},
		{3.4129, 4.9317, 8.7643, 1.3142},
		{1.2371, 4.9998, 10.6721, 0.0147}
	};
	double b[N]={9.5342, 6.3941, 18.4231, 16.9237};
	double m[N][N];
	double x[N];
	int i,j,k,l;
	double max;
	double p;
	for(k=0;k<N-1;k++)
	{
		max=a[k][k];
		l=k;
		for(i=k+1;i<N;i++)//选主元
		{
			if(fabs(a[i][k])>fabs(max))
			{
				max=a[i][k];
				l=i;
			}
		}
		if(fabs(max)*100000<1)
			cout<<"det A=0"<<endl;
		else{
			if(l!=k)//交换行
			{
				p=b[l];
				b[l]=b[k];
				b[k]=p;
				for(j=0;j<N;j++)
				{
					p=a[l][j];
					a[l][j]=a[k][j];
					a[k][j]=p;
				}
				display(a,b);
				cout<<endl;
			}
			for(i=k+1;i<N;i++)//消元计算
			{
				m[i][k]=a[i][k]/a[k][k];
				for(j=k+1;j<N;j++)
				{
					a[i][j]-=m[i][k]*a[k][j];
				}
				b[i]-=m[i][k]*b[k];
			}
			display(a,b);
			cout<<endl;
		}
	}
	if(fabs(a[N-1][N-1])*100000<1)
		cout<<"det A=0"<<endl;
	//回代求解
	else{
		x[N-1]=b[N-1]/a[N-1][N-1];
		for(i=N-2;i>=0;i--)
		{
			x[i]=b[i];
			for(j=i+1;j<N;j++)
			{
				x[i]-=a[i][j]*x[j];
			}
			x[i]/=a[i][i];
		}
	}
	//输出方程组的解
	cout<<"方程的解依次为："<<endl;
	for(i=0;i<N;i++)
	{
		cout<<"x["<<i+1<<"]="<<x[i]<<endl;
	}
	system("pause");
}