学习记录-高斯消元

高斯消元

就是解方程。给n个未知数与n个多元一次方程。即n个方程左边的系数组成一个系数矩阵a[][]，右边的值组成一个矩阵x[]。高斯消元可以在O(N^3)的时间内求出所有未知数的值。具体方法很简单，就是把这n个一次方程用过加减尽可能的消掉未知数，最后求出所有未知数的解。
例如：

那么就有a[][]矩阵和x[]矩阵分别是前三列和最后一列：

**即每一次第（i+1）行的系数减去第（i）行的系数的一定的倍数，来消去一个未知数（即把这个未知数的系数化为0），然后更新a数组第（i+1）行中所有的系数和x数组中的答案，然后从第2行遍历到最后一行，最后的得到的x数组就是答案。**这里还需要注意的是，为了减少误差，每一次的参照行，即目前的第一行，应是第一个系数最大那一行。

模板

其中n与m 表示方程的个数和未知数的个数，即a矩阵的行和列。一般来说n==m

//#include<bits/stdc++.h>
#include<algorithm>
#include<complex>
#include<iostream>
#include<iomanip>
#include<ostream>
#include<cstring>
#include<string.h>
#include<string>
#include<cstdio>
#include<cctype>
#include<vector>
#include<cmath>
#include<queue>
#include<set>
#include<stack>
#include<map>
#include<cstdlib>
#include<time.h>
#include<ctime>
#include<bitset>
#define pb push_back
#define _fileout freopen("out.txt","w",stdout)
#define _filein freopen("in.txt","r",stdin)
#define ok(i) printf("ok%d\n",i)
using namespace std;
typedef double db;
typedef long long ll;
const double PI = acos(-1.0);
const ll MOD=1e9+7;
const ll NEG=1e9+6;
const int MAXN=2e2+10;
const int INF=0x3f3f3f3f;
const ll ll_INF=9223372036854775807;
const double eps=1e-9;
ll qm(ll a,ll b){ll ret = 1;while(b){if (b&1)ret=ret*a%MOD;a=a*a%MOD;b>>=1;}return ret;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return (a*b)/gcd(a,b);}
double a[MAXN][MAXN];//0 base
double x[MAXN];
int n,m;//方程的数量和未知数的个数,相当于m，n。即行和列
int Gauss()
{
    int i,j,k,col,max_r;
    for(k=0,col=0;k<n&&col<m;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<n;i++)//找目前a[k][col]最大的来作为第一个，便于减少误差
            if(fabs(a[i][col])>fabs(a[max_r][col]))
                max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r)//找到之后把它与目前这一行(第k行)交换.
        {
            for(j=col;j<m;j++)
                swap(a[k][j],a[max_r][j]);
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<m;j++)a[k][j]/=a[k][col];//在这里直接把倍数赋给a[k][j]，下面更新时直接乘即可
        a[k][col]=1;
        for(i=0;i<n;i++)
            if(i!=k)
            {
                x[i]-=x[k]*a[i][col];//更新x[i]
                for(j=col+1;j<m;j++)a[i][j]-=a[k][j]*a[i][col];//更新系数
                a[i][col]=0;//这一列的其他的行更新为0
            }
    }
    return 1;//最后返回的数组中，a数字时一个对角线为1，其余位置为0的矩阵，x数组是相对应的答案。
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<m;j++)
        {
            scanf("%lf",&a[i][j]);
        }
        scanf("%lf",&x[i]);
    }
    if(Gauss())
        for(int i=0;i<m;i++)
        {
            printf("%.2f ",x[i]);
        }
    return 0;
}