洛谷 P3389：[模板] 高斯消元法 ← 高斯-约旦消元法（Gauss-Jordan Elimination）

​【题目来源】
https://www.luogu.com.cn/problem/P3389

【题目描述】
给定一个线性方程组，对其求解。

【输入格式】
第一行，一个正整数 n。
第二至 n+1 行，每行 n+1 个整数，为 a1,a2,…,an 和 b，代表一组方程。

【输出格式】
共 n 行，每行一个数，第 i 行为 xi（四舍五入保留 2 位小数）。
如果不存在唯一解或无解，在第一行输出 No Solution。

【输入样例】
3
1 3 4 5
1 4 7 3
9 3 2 2​​​​​​​

【输出样例】
-0.97
5.18
-2.39

【数据范围】
1≤n≤100，|ai|≤10^4，|b|≤10^4。
保证数据若有解则所有解均满足 |xi|≤10^3，且 xi±10^(-6) 和 xi 四舍五入后的结果相同（即不会因为较小的精度误差导致四舍五入后的结果不同）。

【算法分析】
● 高斯消元法（Gaussian Elimination）​ 是一种用于求解线性方程组的经典算法，也是线性代数中最基础的工具之一。它通过一系列初等变换，将增广矩阵化为行阶梯形矩阵或行最简形矩阵，从而求出方程组的解。

● 高斯消元法的本质是系统地消除未知数，将复杂的方程组转化为更容易求解的形式。整个过程分为两个阶段：
（1）前向消元（Forward Elimination）：将矩阵化为上三角矩阵（行阶梯形）。
（2）回代求解（Back Substitution）：从最后一个方程开始，逐个求解未知数。​​​​​​​

● 核心代码解析
（一）选主元

int idx=r;
for(int i=r; i<n; i++) {
    if(fabs(a[i][c])>fabs(a[idx][c])) idx=i;
}

从第 r 行开始，往下找第 c 列上绝对值最大的那个元素。
idx 记录绝对值最大的行号。
为什么要找绝对值最大的？—— 减小浮点误差。

（二）判断是否跳过该列？

if(fabs(a[idx][c])<eps) continue;

如果找到的最大值都接近于 0，说明这一列全是 0。
这个未知数在当前剩余的方程中不存在，直接跳过该列（它可能是自由变量）。
continue 跳过后面的代码，进入下一列，但 r 不变（还没找到新的主元）。

（三）交换到当前行

swap(a[r],a[idx]);

把找到的主元行（第 idx 行）换到第 r 行。
swap(a[r], a[idx]) 交换的是整行，包括常数项。

（四）归一化​​​​​​​

double div=a[r][c];
for(int j=c; j<=n; j++) a[r][j]/=div;

div 就是主元的值。
把第 r 行的所有元素（从第 c 列到第 n 列，含常数项）都除以 div。
结果：主元变成了 1，这一行其他元素相应缩放。

（五）消去其他行​​​​​​​

for(int i=0; i<n; i++) {
    if(i!=r && fabs(a[i][c])>eps) {
        double val=a[i][c];
        for(int j=c; j<=n; j++) {
            a[i][j]-=val*a[r][j];
        }
    }
}

对所有行（不只是下面的行），如果第 c 列不为 0，就消掉它。
val = a[i][c] 是第 i 行第 c 列的当前值。
a[i][j] -= val * a[r][j] 就是用第 r 行（归一化后主元为 1）去消。
消完后，第 i 行的第 c 列变成 0。
注意：这里消的是所有行（包括上面的），这是高斯-约旦消元法的特点，直接得到对角矩阵，不用回代。

【算法代码】

#include <bits/stdc++.h>
using namespace std;

const int N=105;
const double eps=1e-8;
double a[N][N];
int n;

int gauss() {
    int r=0;
    for(int c=0; c<n; c++) {
        //1.Select pivot
        int idx=r;
        for(int i=r; i<n; i++) {
            if(fabs(a[i][c])>fabs(a[idx][c])) idx=i;
        }
        if(fabs(a[idx][c])<eps) continue;
        swap(a[r],a[idx]);

        //2.Normalize the current row, and the pivot becomes 1
        double div=a[r][c];
        for(int j=c; j<=n; j++) a[r][j]/=div;

        //3.Eliminate all rows below
        for(int i=0; i<n; i++) {
            if(i!=r && fabs(a[i][c])>eps) {
                double val=a[i][c];
                for(int j=c; j<=n; j++) {
                    a[i][j]-=val*a[r][j];
                }
            }
        }
        r++;
    }


    for(int i=r; i<n; i++) { //No solution
        if(fabs(a[i][n])>eps) return -1;
    }
    if(r<n) return 0; //Infinite solutions
    return 1; //Unique solution
}

int main() {
    cin>>n;
    for(int i=0; i<n; i++) {
        for(int j=0; j<=n; j++) {
            cin>>a[i][j];
        }
    }

    int t=gauss();
    if(t==-1 || t==0) cout<<"No Solution\n";
    else {
        for(int i=0; i<n; i++) {
            printf("%.2lf\n",a[i][n]);
        }
    }
    return 0;
}

/*
in:
3
1 3 4 5
1 4 7 3
9 3 2 2

out:
-0.97
5.18
-2.39
*/

【参考文献】
https://www.luogu.com.cn/problem/solution/P3389
 

​