[板子] 高斯消元法

Gauss Elimination

流程

(i) 把系数矩阵化成上三角矩阵
    1. 对每一行,选出还未处理过的行中不为0的首项的绝对值最大的一行,和当前行交换。
    2. 把当前行的每一项除以不为0的首项。
    3. 把未处理过的每一行的首项通过使此行减当前行的某倍消去。
(ii) 从最后一行把已知的解带回上一行的方程,得出上一行的解。

Eg. 矩阵

解过程如下:

接下来回带求解:
显然$z = 3$,带入上式得

得$y = 2$。最后把y和z的值带入第一行:

得$x = 1$。


Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <cstdio>
#include <algorithm>
#include <cmath>
const int N = 150,inf = 1e9;
struct Gauss
{
int n;double a[N][N],ret[N];
bool gauss()
{
for(int i = 1;i <= n;i++)
{
int m = i;
for(int j = i;j <= n;j++) if(std::fabs(a[m][i]) < std::fabs(a[j][i])) m = j;
for(int j = 1;j <= n + 1;j++) std::swap(a[i][j],a[m][j]);
for(int j = n + 1;j;j--) a[i][j] /= a[i][i];
for(int j = i + 1;j <= n;j++)
{
double l = a[j][i] / a[i][i];
for(int k = i;k <= n + 1;k++) a[j][k] -= a[i][k] * l;
}
}
for(int i = n;i;i--)
{
for(int j = n;j >= i + 1;j--) a[i][n + 1] -= ret[j] * a[i][j];
ret[i] = a[i][n + 1];
}
for(int i = 1;i <= n;i++) if(ret[i] != ret[i]) return 0;
return 1;
}
}G;
int main()
{
int n;scanf("%d",&n);G.n = n;
for(int i = 1;i <= n;i++) for(int j = 1;j <= n + 1;j++) {int x;scanf("%d",&x);G.a[i][j] = 1.0 * x;}
if(!G.gauss()) printf("No Solution\n");
else for(int i = 1;i <= n;i++) printf("%.2f\n",G.ret[i]);
return 0;
}