浅谈高斯消元
模板题：P3389 【模板】高斯消元法
本文采用高斯-约旦消元法（英语：Gauss-Jordan Elimination）。
相较传统的高斯消元，它不需要回代的过程，所以比较方便。
高斯消元通常用于求解 \(n\) 元线性方程组的解（或返回无解），下面是 \(n=3\) 的例子 \[ \begin{cases} 3 x+2 y+z=10 \\[6pt]5 x+y+6z=25 \\[6pt]2 x+3 y+4 z=20 \end{cases} \] 我们可以把它写成矩阵的形式 \[ \begin{bmatrix} &3&2&1&10& \\[6pt]&5&1&6&25& \\[6pt]&2&3&4&20& \end{bmatrix} \] 一般地，考虑 \(n\) 元线性方程组 \[ \begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\[6pt]a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\[6pt]\quad\vdots \\[6pt]a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_2 \end{cases} \] 其矩阵形式为 \[ \begin{bmatrix} &a_{11}&a_{12}&\cdots&a_{1n}&b_1& \\[6pt]&a_{21}&a_{22}&\cdots&a_{2n}&b_2& \\[6pt]&\vdots&\vdots&\ddots&\vdots&\vdots \\[6pt]&a_{n1}&a_{n2}&\cdots&a_{nn}&b_n& \end{bmatrix} \] 依次考虑每个 \(a_{ii}\) ，将除第 \(i\) 行的每一行第 \(i\) 列都变为 \(0\) 。
具体操作为，对于每个 \(a_{jk}\) ，使其变为 \(a_{jk} - a_{ik} \times \frac{a_{ji}}{a_{ii}}\) 。
最后 \(x_i\) 的解就是 \(\frac{b_i}{a_{ii}}\) ，且最后的矩阵应该长这样（空白处为 \(0\) ） \[ \begin{bmatrix} &a_{11}&&&&b_1& \\[6pt]&&a_{22}&&&b_2& \\[6pt]&&&\ddots&&\vdots \\[6pt]&&&&a_{nn}&b_n& \end{bmatrix} \]
实现时有一些细节，每次把第 \(i \sim n\) 行的第 \(i\) 列的系数最大值找出来
把它与第 \(i\) 行交换，这样可以保证算法的稳定性。（如果最后 \(a_{ii}=0\) 则无解）
时间复杂度 \(\mathcal{O}(n^3)\)
代码：
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(150))
const double eps = 1e-10;
int dcmp(int x)
{
if(abs(x) <= eps) return 0;
return x > eps ? 1 : -1;
}
int n; double a[N][N];
signed main()
{
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
// freopen("check.in","r",stdin);
// freopen("check.out","w",stdout);
cout << fixed << setprecision(2);
cin >> n;
for(int i=1; i<=n; i++)
for(int j=1; j<=n+1; j++)
cin >> a[i][j];
for(int i=1; i<=n; i++) // col
{
int p = i;
for(int j=i+1; j<=n; j++)
if(abs(a[j][i]) > abs(a[p][i])) p=j; // row
for(int j=1; j<=n+1; j++)
swap(a[i][j],a[p][j]); // swap p to i for proc
if(!dcmp(a[i][i])) return cout << "No Solution\n",0;
for(int j=1; j<=n; j++)
if(j != i)
{
double tmp = a[j][i] / a[i][i];
for(int k=1; k<=n+1; k++)
a[j][k] -= a[i][k] * tmp;
}
}
for(int i=1; i<=n; i++)
cout << a[i][n+1] / a[i][i] << '\n';
return 0;
}参考文献：
