【学习笔记】高斯消元法

P3389 【模板】高斯消元法

数论小白都能看懂的线性方程组及其解法

约旦消元法：

• 首先找到当前系数（i）的最大（主元）的柿子（为了尽可能减少爆精度），并调整到第i行
• 若当前最大项系数为0，No Solution
• 接下来进行对每个系数的非主元行的消元。形成俩全0三角形

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define dwn(i,a,b) for(int i=a;i>=b;--i)
template <typename T> inline void rd(T &x){x=0;char c=getchar();int f=0;while(!isdigit(c)){f|=c=='-';c=getchar();}while(isdigit(c)){x=x*10+(c^48);c=getchar();}x=f?-x:x;}

int n;
double a[110][110];

int main(){
    #ifdef WIN32
    freopen("xiaoyuan.txt","r",stdin);
    #endif
    rd(n);
    rep(i,1,n)
        rep(j,1,n+1)
            scanf("%lf",&a[i][j]);
    rep(i,1,n){//枚举每一列的系数
        int max=i;//假定这一列最大的系数是第i行的；
        rep(j,i+1,n)
            if(fabs(a[j][i])>fabs(a[max][i]))
                max=j;
        rep(k,1,n+1)
            swap(a[max][k],a[i][k]);//这样每次第i个未知数最大的系数都在第i行，可以保证俩个0三角形
        if(!a[i][i]){//最大值等于0则说明该列都为0，肯定无解
            printf("No Solution");
            return 0;
        }
        rep(j,1,n){//每一行都减去一个数
            if(i!=j){//不是主元那一行
                double tmp=a[j][i]/a[i][i];//注意这里是浮点数！！！！
                rep(k,i,n+1)
                a[j][k]-=a[i][k]*tmp;
            }
        }
    }
    rep(i,1,n)
        printf("%.2lf\n",a[i][n+1]/a[i][i]);
    return 0;
}
/*
3
1 3 4 5
1 4 7 3
9 3 2 2
*/
/*
-0.97
5.18
-2.39
*/