【高斯消元】浮点高斯消元

浮点数高斯消元
【问题描述】
给出一个线性方程组，有n个方程组，m个未知数。解这个线性方程组。

【输入格式】
第1行：2个整数n和m，(n, m <=400，且n不一定等于m)
接下来n行，每行m+1个整数，表示一个方程的m个未知数的系数和常数

【输出格式】
如果无解，输出“No solution”。
如果有唯一解，输出m行，每行一个未知数的值，保留到小数点第4位。格式见样例。
如果有无穷解，输出m行，如果未知数有确定解，直接输出。如果是自由变元，输出“xx isfree number”

【输入样例1】
3 3
2 -1 3 1
4 2 5 4
2 0 2 6
【输出样例1】
X[1] = 9.0000
X[2] = -1.0000
X[3] = -6.0000

【输入样例2】
3 3
2 -1 3 1
4 -2 5 4
2 -1 4 -1

【输出样例2】
X[1] not determined
X[2] not determined
X[3] = -2.0000

【输入样例3】
3 4
5 -1 2 1 7
2 1 4 -2 1
1 -3 -6 5 0

【输出样例3】
No solution

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define SF scanf
#define PF printf
#define MAXN 410
#define EPS 1e-8
using namespace std;
double a[MAXN][MAXN],x[MAXN];
int n,m,rank,f[MAXN];
void init(){
    memset(a,0,sizeof a);
    for(int i=0;i<n;i++)
        for(int j=0;j<=m;j++)
            SF("%lf",&a[i][j]);
}
void gauss(){
    int r,c,maxr;
    for(r=0,c=0;r<n&&c<m;r++,c++){
        maxr=r;
        for(int i=r+1;i<n;i++)
            if(fabs(a[i][c])>fabs(a[maxr][c]))
                maxr=i;
        if(fabs(a[maxr][c])<EPS){
            r--;
            continue;
        }
        if(maxr!=r)
            for(int i=c;i<=m;i++)
                swap(a[r][i],a[maxr][i]);
        for(int i=0;i<n;i++)
            if(i!=r&&fabs(a[i][c])>EPS)
                for(int j=m;j>=c;j--)
                    a[i][j]-=a[r][j]/a[r][c]*a[i][c];
    }
    rank=r;
}
bool check(){
    int cnt,pos;
    for(int i=rank;i<n;i++)
        if(fabs(a[i][m])>EPS)
            return 0;
    for(int i=0;i<m;i++)
        f[i]=1;
    for(int i=rank-1;i>=0;i--){
        cnt=0;
        for(int j=0;j<m;j++)
            if(fabs(a[i][j])>EPS&&f[j]){
                cnt++;
                pos=j;
            }
        if(cnt==1){
            f[pos]=0;
            x[pos]=a[i][m]/a[i][pos];
        }
    }
    return 1;
}
void print(){
    int i;
    if(rank<m){
        PF("Multiple solution, free_num: %d\n",m-rank);
        for(int i=0;i<m;i++){
            if(f[i])
                PF("X[%d] not determined\n",i+1);
            else
                PF("X[%d] = %.4f\n",i+1,x[i]);
        }
    }
    else{
        for(int i=0;i<n;i++)
            PF("X[%d] = %.4f\n",i+1,x[i]);
    }
}
int main(){
    while(SF("%d%d",&n,&m)!=EOF){
        init();
        gauss();
        if(check())
            print();
        else
            PF("No solution\n");
    }
}