【模板】高斯消元法

题目背景

Gauss 消元

题目描述

给定一个线性方程组，对其求解

输入格式

第一行，一个正整数 $n$

第二至 $n+1$ 行，每行 $n+1$ 个整数，为$ a_1, a_2 \cdots a_n$ 和 $b$，代表一组方程。

输出格式

共 $n$ 行，每行一个数，第 $i$ 行为 $x_i$ （保留 2 位小数）

如果不存在唯一解，在第一行输出 No Solution.

样例 #1

样例输入 #1

样例输出 #1

1
2
3
-0.97
5.18
-2.39

提示

$1 \leq n \leq 100, \left | a_i \right| \leq {10}^4 , \left |b \right| \leq {10}^4 $

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
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <bits/stdc++.h>
using namespace std;

int sz;
bool no_solution = false;
const double eps = 1e-9;
bool IsZero(double v) { return abs(v) < eps; }
enum GAUSS_MODE { DEGREE, ABS };
template <typename T>
void GaussianElimination(vector<vector<T>>& a, int limit,
                         GAUSS_MODE mode = DEGREE) {
    if (a.empty() || a[0].empty()) return;
    int h = static_cast<int>(a.size());
    int w = static_cast<int>(a[0].size());
    vector<int> deg(h);
    for (int i = 0; i < h; i++)
        for (int j = 0; j < w; j++) deg[i] += !IsZero(a[i][j]);
    int r = 0;
    for (int c = 0; c < limit; c++) {
        int id = -1;
        for (int i = r; i < h; i++)
            if (!IsZero(a[i][c]) &&
                (id == -1 || (mode == DEGREE && deg[i] < deg[id]) ||
                 (mode == ABS && abs(a[id][c]) < abs(a[i][c]))))
                id = i;
        if (id == -1) continue;
        if (id > r) {
            swap(a[r], a[id]);
            swap(deg[r], deg[id]);
            for (int j = c; j < w; j++) a[id][j] = -a[id][j];
        }
        vector<int> nonzero;
        for (int j = c; j < w; j++)
            if (!IsZero(a[r][j])) nonzero.push_back(j);
        T inv_a = 1 / a[r][c];
        for (int i = r + 1; i < h; i++) {
            if (IsZero(a[i][c])) continue;
            T coeff = -a[i][c] * inv_a;
            for (int j : nonzero) {
                if (!IsZero(a[i][j])) deg[i]--;
                a[i][j] += coeff * a[r][j];
                if (!IsZero(a[i][j])) deg[i]++;
            }
        }
        ++r;
    }
    for (r = h - 1; r >= 0; r--) {
        for (int c = 0; c < limit; c++) {
            if (!IsZero(a[r][c])) {
                T inv_a = 1 / a[r][c];
                for (int i = r - 1; i >= 0; i--) {
                    if (IsZero(a[i][c])) {
                        continue;
                    }
                    T coeff = -a[i][c] * inv_a;
                    for (int j = c; j < w; j++) {
                        a[i][j] += coeff * a[r][j];
                    }
                }
                break;
            }
        }
    }
}

template <typename T>
vector<T> SolveLinearSystem(vector<vector<T>> a, const vector<T>& b, int w) {
    no_solution = false;
    int h = static_cast<int>(a.size());
    assert(h == static_cast<int>(b.size()));
    if (h > 0) {
        assert(w == static_cast<int>(a[0].size()));
    }
    for (int i = 0; i < h; i++) {
        a[i].push_back(b[i]);
    }
    GaussianElimination(a, w);
    vector<T> x(w, 0);
    for (int i = 0; i < h; i++) {
        for (int j = 0; j < w; j++) {
            if (!IsZero(a[i][j])) {
                x[j] = a[i][w] / a[i][j];
                break;
            }
            if (j == w - 1) {
                no_solution = true;
                return vector<T>();
            }
        }
    }
    return x;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("data/in.txt", "r", stdin);
    freopen("data/out.txt", "w", stdout);
#endif
    ios::sync_with_stdio(false);
    int n;
    double x;
    cin >> n;
    vector<vector<double>> a(n);
    vector<double> b;
    for (int i = 0; i < n; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> x;
            a[i].push_back(x);
        }
        cin >> x;
        b.push_back(x);
    }
    vector<double> c = SolveLinearSystem(a, b, n);
    if (no_solution)
        cout << "No Solution";
    else
        for (int i = 0; i < c.size(); i++) printf("%.2lf\n", c[i]);
    return 0;
}