首页 > 技术文章 > 高斯消元解线性方程组

IamIron-Man 2020-09-24 11:55 原文

线性方程组形如,高斯消元可以用O(n^3)的复杂度解出x={}。
二维vector数组定义后,只能通过push_back一个个读入。
C++代码:

#include<iostream>
#include<iomanip>
#include<cstdlib>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<queue>
#include<vector>
#include<cmath> 
#include<map> 
#define ll long long
#define pii pair<int,int>
#define mpr make_pair
#define inf 0x3f3f3f3f 
#define INF 0x3f3f3f3f3f3f3f3f
#define clr(x) memset(x,0,sizeof x)
#define rep2(i, x) for(register int i = h[x]; i; i = e[i].nxt)
#define rep(i,x,y) if((x)<=(y))for(register int i=(x);i<=(y);i++)
#define rrep(i,x,y) if((x)>=(y))for(register int i=(x);i>=(y);i--)
#define endl "\n" 
#define TT() int t;cin>>t;while(t--) 
#define PI acos(-1)
#define lson l , m , rt << 1
#define rson m + 1 , r , rt << 1 | 1
const double EPS=1e-8;
using namespace std;
typedef vector<double> vec;
typedef vector<vec> mat;

vec gauss_jordan(const mat &A,const vec &b) {
	int n=A.size();
	mat B(n,vec(n+1));
	rep(i,0,n-1)
	  rep(j,0,n-1)
	    B[i][j]=A[i][j];
	rep(i,0,n-1)
		B[i][n]=b[i];
	rep(i,0,n-1)
	{
		int pivot=i;
		rep(j,i,n-1)
			if (abs(B[j][i])>abs(B[pivot][i])) 
				pivot=j;
		swap(B[i],B[pivot]);
		if (abs(B[i][i])<EPS) return vec();
		rep(j,i+1,n)
		  B[i][j]/=B[i][i];
		rep(j,0,n-1)
		{
		  if (i!=j)
		  {
		  	rep(k,i+1,n)
		  	  B[j][k]-=B[j][i]*B[i][k];
		  }
		}
	}
	vec x(n);
	rep(i,0,n-1)
	  x[i]=B[i][n];
	return x;
}

int main(){
	std::ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	mat a;
	int n,m,t2;
	vec b,x;
	cin>>n>>m;
	rep(i,0,n-1)
	{
		vec t1;
	    rep(j,0,m-1)
	    {
	  	  cin>>t2;
	  	  t1.push_back(t2);
	    }
	    a.push_back(t1);
	}
	rep(i,0,n-1)
	{
	  cin>>t2;
	  b.push_back(t2);
    }
	x=gauss_jordan(a,b);
	rep(i,0,x.size()-1)
	  cout<<x[i]<<' ';
	return 0;
}

推荐阅读