高斯消元法

高斯消元法

#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>

using namespace std;

const double EPS=1e-8;

typedef vector<double>vec;
typedef vector<vec>mat;


vec Gauss(const mat &A,const vec &b)
{
	int n=A.size();//n个未知数
	mat B(n,vec(n+1));//初始化增广矩阵
	for(int i=0;i<n;i++)
	  for(int j=0;j<n;j++)
	    B[i][j]=A[i][j];
	
	for(int i=0;i<n;i++)//常数矩阵 
	  B[i][n]=b[i];
	
	//把正在处理的未知数系数的绝对值最大的式子换到第i行
	for(int i=0;i<n;i++)
	{
		int pos=i;
		for(int j=i;j<n;j++)
		{
			if(abs(B[j][i])>abs(B[pos][i]))
			  pos=j;
		}
		swap(B[i],B[pos]);
		
		//无解或者有无穷解
		if(abs(B[i][i])<EPS) return vec();
		
		//把正在处理的未知数的系数变为1
		double tmp=B[i][i];
		for(int j=i;j<=n;j++) B[i][j]/=tmp;
		
		for(int j=0;j<n;j++)
		{
			if(i!=j)
			{
				for(int k=i+1;k<=n;k++)
				  B[j][k]-=B[j][i]*B[i][k];
			}
		 } 
	}
	vec x(n);
	for(int i=0;i<n;i++) x[i]=B[i][n];
	return x;
}

int main()
{
	int n;
	scanf("%d",&n);
	mat A(n,vec(n));
	vec b(n);
	//输入系数矩阵 
	for(int i=0;i<n;i++)
	  for(int j=0;j<n;j++)
	    scanf("%lf",&A[i][j]);
	
	//输入常数矩阵
	for(int i=0;i<n;i++)
	  scanf("%lf",&b[i]);
	
	vec res=Gauss(A,b);
	for(int i=0;i<res.size();i++)
	  printf("%lf ",res[i]);
	return 0; 
}