高斯消元解线性方程组

    在线性代数中我们会经常遇到诸如求以下解线性方程组的问题：

x1*a11+x2*a12+...xn*a1n=b1

x1*a21+x2*a22+...xn*a2n=b2

...

x1*am1+x2*am2+...xn*amn=bm

    我们可以用高斯消元法来求解上述方程组的解，具体实现代码如下：

整型系数代码如下：

/*
Author:Ibsen
Date:2015.11.16
*/
#include <cstdio>
#include <iostream>
using namespace std;
const int M=100; //处理的矩阵阶数的最大值
int a[M][M]; //系数矩阵和增广矩阵
int x[M]; //存放最后的解
int equ,var; //equ是系数阵的行数,var是系数阵的列数
void Debug()
{
	cout<<"消元结果为："<<endl;
	for(int i=0;i<equ;i++)
	{
		for(int j=0;j<var+1;j++)
			cout<<a[i][j]<<" ";
		cout<<endl;
	}
}
int abs(int x)
{
    if(x>=0) return x;
    else return -x;
}
inline int Gcd(int a,int b)
{
	return !b?a:Gcd(b,a%b);
}
inline int Lcm(int a,int b)
{
	return a*b/Gcd(a,b);
}
void swap(int &a,int &b)
{
	int tmp=a;
	a=b;
	b=tmp;
}
int Gauss()
{
	int k,col=0; //k表示当前处理的行,col表示当前处理的列
	for(k=0;k<equ&&col<var;k++,col++)
	{
		int max_r=k;
		for(int i=k+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[max_r][col])) //这里要找绝对值最大的列，避免代码尾部样例第二组（只有负数和0，最大值为0）的情况的情况
			  max_r=i;
		if(max_r!=k)
		{
			for(int i=k;i<var+1;i++)
				swap(a[k][i],a[max_r][i]);
		}
		if(!a[k][col]) //(第k行往后)第col列元素全为0,直接进入下一列
		{
			k--;
			continue;
		}
		for(int i=k+1;i<equ;i++)
		{
			if(a[i][col])
			{
				int lcm=Lcm(a[i][col],a[k][col]);
				int ta=lcm/a[i][col], tb=lcm/a[k][col];
				for(int j=col;j<var+1;j++)
					a[i][j]=a[i][j]*ta-a[k][j]*tb;
			}
		}
	}
	for(int i=k;i<equ;i++)
		if(a[i][col]) return -1; //无解
	Debug();
	for(int i=0;i<equ;i++)
	{
		if(!a[i][i])
		{
			int j=i+1;
			for(;j<var;j++)
				if(a[i][j]) break;
			if(j==var) break;
			for(int k=0;k<equ;k++)
				swap(a[k][i],a[k][j]);
		}
	}
	Debug();
	cout<<"col="<<col<<endl;
	if(var-k==0) //有唯一解：回带求解
	{
		for(int i=k-1;i>=0;i--) //找出最后一个非零行
		{
			int tmp=a[i][var];
			for(int j=i+1;j<var;j++)
				if(a[i][j])
				{
					tmp=tmp-a[i][j]*x[j];
				}
			x[i]=tmp/a[i][i];
		}
	}
}
void display()
{
	cout<<"解向量X为：";
	for(int i=0;i<var;i++)
		cout<<x[i]<<" ";
	cout<<endl;
}
int main()
{
	while(cin>>equ>>var)
	{
		for(int i=0;i<equ;i++)
			for(int j=0;j<var+1;j++)
				cin>>a[i][j];
		if(Gauss()==-1) puts("No solution!");
		else display();
	}
	return 0;
}
/*
2 2
4 5 13
-6 3 -9
解向量为(2,1)
4 4
2 -1 0 0 1
-1 2 -1 0 0
0 -1 2 -1 0
0 0 -1 2 6
解向量为(2,3,4,5)
*/

高精度实现代码如下：

/*
Author:Ibsen
Date:2015.11.16
*/
#include <cstdio>
#include <iostream>
using namespace std;
const int M=100; //处理的矩阵阶数的最大值
const double esp=1e-8; //精度
double a[M][M]; //系数矩阵和增广矩阵
double x[M]; //存放最后的解
int equ,var; //equ是系数阵的行数,var是系数阵的列数
double abs(double x)
{
	if(x>0) return x;
	else return -x;
}
void Debug()
{
	cout<<"消元结果为："<<endl;
	for(int i=0;i<equ;i++)
	{
		for(int j=0;j<var+1;j++)
			cout<<a[i][j]<<" ";
		cout<<endl;
	}
}
void swap(double &a,double &b)
{
	double tmp=a;
	a=b;
	b=tmp;
}
int Gauss()
{
	int k,col=0; //k表示当前处理的行,col表示当前处理的列
	for(k=0;k<equ&&col<var;k++,col++)
	{
		int max_r=k;
		for(int i=k+1;i<equ;i++)
			if(abs(a[i][col])>abs(a[max_r][col])) max_r=i;
		if(max_r!=k)
		{
			for(int i=k;i<var+1;i++)
				swap(a[k][i],a[max_r][i]);
		}
		if(abs(a[k][col])<esp) //(第k行往后)第col列元素全为0,直接进入下一列
		{
			k--;
			continue;
		}
		for(int i=k+1;i<equ;i++)
		{
			if(abs(a[i][col])>esp)
			{
				double tmp=a[k][col]/a[i][col];
				for(int j=col;j<var+1;j++)
					a[i][j]=a[i][j]*tmp-a[k][j];
			}
		}
	}
	for(int i=k;i<equ;i++)
		if(abs(a[i][col])<esp) return -1; //无解
	//Debug();
	for(int i=0;i<equ;i++)
	{
		if(abs(a[i][i])<esp)
		{
			int j=i+1;
			for(;j<var;j++)
				if(abs(a[i][j])>esp) break;
			if(j==var) break;
			for(int k=0;k<equ;k++)
				swap(a[k][i],a[k][j]);
		}
	}
	Debug();
	cout<<"col="<<col<<endl;
	if(var-k==0) //有唯一解：回带求解
	{
		for(int i=k-1;i>=0;i--) //找出最后一个非零行
		{
			double tmp=a[i][var];
			for(int j=i+1;j<var;j++)
				if(abs(a[i][j])>esp)
				{
					tmp=tmp-a[i][j]*x[j];
				}
			x[i]=tmp/a[i][i];
		}
	}
}
void display()
{
	cout<<"解向量X为：";
	for(int i=0;i<var;i++)
		cout<<x[i]<<" ";
	cout<<endl;
}
int main()
{
	while(cin>>equ>>var)
	{
		for(int i=0;i<equ;i++)
			for(int j=0;j<var+1;j++)
				cin>>a[i][j];
		if(Gauss()==-1) puts("No solution!");
		else display();
	}
	return 0;
}
/*
3 3
1.0 2.0 3.0 6.6
2.0 3.0 4.0 9.9
1.0 3.0 2.0 6.6
解向量为(1.1,1.1,1.1)
*/