A = LU 分解
$L$ 矩阵是一个下三角矩阵，$U$ 矩阵是一个上三角矩阵。这种表示方式相对于 $EA = U$ 而言更加合适。
$$
\quad
$$
高斯消元的时间复杂度是 $O(n ^ 3)$ 级别的，编程难度不算很大。
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
39
40
41
| #include <bits/stdc++.h>
const int N = 150;
const double eps = 1e-7;
using namespace std;
int n;
double ans[N], a[N][N];
int main() {
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
freopen("in","r", stdin);
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++) {
int pos = i;
for (int k = i + 1; k <= n; k++) {
if (fabs(a[i][i]) < fabs(a[k][i])) pos = k;
}
if (fabs(a[pos][i]) < eps) {
cout << "No Solution" << endl;
return 0;
}
if (pos != i) swap(a[i], a[pos]);
double div = a[i][i];
for (int j = i; j <= n + 1; j++) a[i][j] /= div;
for (int k = i + 1; k <= n; k++) {
div = a[k][i];
for (int l = i; l <= n + 1; l++) a[k][l] -= div * a[i][l];
}
}
ans[n] = a[n][n + 1];
for (int i = n - 1; i >= 1; i--) {
ans[i] = a[i][n + 1];
for (int j = i + 1; j <= n; j++) ans[i] -= ans[j] * a[i][j];
}
for (int i = 1; i <= n; i++) cout << fixed << setprecision(2) << ans[i] << endl;
return 0;
}
|
对于消元过程中碰到 pivot 为 0 的情况，需要进行行交换。行交换等价于在原矩阵的左边乘一个置换矩阵（permutation matrix）。因此，$A = LU$ 分解写成 $PA = LU$ 的形式将更完整。
置换矩阵 $P$ 的一个重要性质： $P ^ {-1} = P ^ T$。