【数学】高斯消元法

列主元GaussJordan消元法

struct Matrix {

    static const int MAXN = 505;
    static const int MAXM = 505;

    int n, m;
    long double A[MAXN][MAXM];
    long double x[MAXM];
    int pos[MAXM];

    Matrix() {
        ms(A);
    }

    void show() {
        printf("Matrix n=%d m=%d
", n, m);
        for (int i = 1; i <= n; ++i) {
            for (int j = 1; j <= m; ++j)
                printf("%6.2f ", A[i][j]);
            printf("| %6.2f
", A[i][m + 1]);
        }
    }

    // 列主元GaussJordan消元法
    int GaussJordan() {
        // 把矩阵化简为行最简形矩阵
        int row = 0; // 当前已有的非零行的数量
        for(int col = 1; col <= m && row <= n; ++col) {
            // 在剩余行中，寻找列主元所在行
            int pivot = row + 1;
            for(int i = pivot + 1; i <= n; ++i) {
                if(abs(A[i][col]) > abs(A[pivot][col]))
                    pivot = i;
            }
            if(abs(A[pivot][col]) <= EPS) {
                pos[col] = -1;  // 在剩余行中，列主元为0，该变量是自由变量
                continue;
            }
            pos[col] = ++row;
            if(pivot != row) {
                // 行交换，行列式的值变成相反数
                for (int j = m + 1; j >= col; --j)
                    swap(A[pivot][j], A[row][j]);
            }

            // 列主元置1
            for (int j = m + 1; j >= col; --j)
                A[row][j] /= A[row][col];

            // 在所有其他行中，消去该变量
            for(int i = 1; i <= n; ++i) {
                if(i == row || abs(A[i][col]) <= EPS)
                    continue;  // 跳过当前行和不存在该变量的行
                for (int j = m + 1; j >= col; --j)
                    A[i][j] -= A[row][j] * A[i][col];
            }
        }
        // 检查无解的情形
        for(int i = 1; i <= n; ++i) {
            if(abs(A[i][m + 1]) <= EPS)
                continue;  // 空行/非空行 = 0值，均不影响
            bool ok = 0;
            for(int j = 1; j <= m; ++j) {
                if(abs(A[i][j]) <= EPS)
                    continue;
                ok = 1;
                break;
            }
            if(!ok)
                return -1;  // 失败：空行 = 非0值，返回-1
        }
        // 回代
        for(int j = 1; j <= m; ++j) {
            if(pos[j] == -1) {
                x[j] = 0; // 自由变量，可以自由取值，默认置0
                continue;
            }
            x[j] = A[pos[j]][m + 1];  // 列主元置1，自由变量均置0的情形
        }
        return row; //成功：返回矩阵的秩
    }

} mat;

Gauss转换为行最简形，传入原本增广矩阵的左侧A（记得m是A的列数，不是增广矩阵的列数），传出一个位置向量，表示原本矩阵的新矩阵的第i行应原本矩阵的第pos[i]行，传出一个系数向量，表示在交换位置之后，位于这个新位置的行所乘的k。
Solve传入一个位置向量，和原本增广矩阵的右侧b，对这个增广矩阵进行求解，解出所有约束变量的取值（自由变量是可以任意取值的）。

【小贴士】

• 未受到方程组约束的自由变量，可以任意取值。对于自由变量的任意一组取值，都可以对约束变量进行适当赋值得到一组解。常见于统计模2意义下的高斯消元的解的个数，注意首先要保证方程组有解。
• 高斯消元法分为两个步骤，首先先把 (Ax=b) 的 (A) 转换成行最简形，复杂度为 (O(n^3)) ，然后代入正确位置的 (b) 进行求解，复杂度为 (O(n)) ，在 CF113D 中， (A) 是不变的，可以利用同一个 (A) 对不同的 (b) 分别解出 (x) 。（也可以使用LU分解法）

两个步骤高斯消元法：

struct Matrix {

    static const int MAXN = 505;
    static const int MAXM = 505;

    int n, m;
    long double A[MAXN][MAXM];
    long double x[MAXM];
    int pos[MAXM];

    Matrix() {
        ms(A);
    }

    void show() {
        printf("Matrix n=%d m=%d
", n, m);
        for (int i = 1; i <= n; ++i) {
            for (int j = 1; j <= m; ++j)
                printf("%6.2f ", (double)A[i][j]);
            printf("| %6.2f
", (double)A[i][m + 1]);
        }
    }

    struct Operation {
        int type;
        int x, y;
        long double k;
    } op;
    vector<Operation> vec;

    // 列主元GaussJordan消元法
    int GaussJordan1() {
        vec.clear();
        // 把矩阵化简为行最简形矩阵
        int row = 0; // 当前已有的非零行的数量
        for(int col = 1; col <= m && row <= n; ++col) {
            // 在剩余行中，寻找列主元所在行
            int pivot = row + 1;
            for(int i = pivot + 1; i <= n; ++i) {
                if(abs(A[i][col]) > abs(A[pivot][col]))
                    pivot = i;
            }
            if(abs(A[pivot][col]) <= EPS) {
                pos[col] = -1;  // 在剩余行中，列主元为0，该变量是自由变量
                continue;
            }
            pos[col] = ++row;
            if(pivot != row) {
                vec.pb({1, pivot, row, 0});
                // 行交换，行列式的值变成相反数
                for (int j = m + 1; j >= col; --j)
                    swap(A[pivot][j], A[row][j]);
            }

            vec.pb({2, row, 0, A[row][col]});
            // 列主元置1
            for (int j = m + 1; j >= col; --j)
                A[row][j] /= A[row][col];

            // 在所有其他行中，消去该变量
            for(int i = 1; i <= n; ++i) {
                if(i == row || abs(A[i][col]) <= EPS)
                    continue;  // 跳过当前行和不存在该变量的行
                vec.pb({3, i, row, A[i][col]});
                for (int j = m + 1; j >= col; --j)
                    A[i][j] -= A[row][j] * A[i][col];
            }
        }
        return row;  //成功：返回矩阵的秩
    }

    int GaussJordan2() {
        for(const Operation &op : vec) {
            int type = op.type, x = op.x, y = op.y;
            long double k = op.k;
            if(type == 1) {
                swap(A[x][m + 1], A[y][m + 1]);
                continue;
            }
            if(type == 2) {
                A[x][m + 1] /= k;
                continue;
            }
            if(type == 3) {
                A[x][m + 1] -= A[y][m + 1] * k;
                continue;
            }
        }
        // 检查无解的情形
        for(int i = 1; i <= n; ++i) {
            if(abs(A[i][m + 1]) <= EPS)
                continue;  // 空行/非空行 = 0值，均不影响
            bool ok = 0;
            for(int j = 1; j <= m; ++j) {
                if(abs(A[i][j]) <= EPS)
                    continue;
                ok = 1;
                break;
            }
            if(!ok)
                return -1;  // 失败：空行 = 非0值，返回-1
        }
        // 回代
        for(int j = 1; j <= m; ++j) {
            if(pos[j] == -1) {
                x[j] = 0; // 自由变量，可以自由取值，默认置0
                continue;
            }
            x[j] = A[pos[j]][m + 1];  // 列主元置1，自由变量均置0的情形
        }
        return 1;  //成功
    }

} mat;

求解行列式：

const double EPS = 1E-9;
int n;
vector<vector<double> > a(n, vector<double>(n));

double det = 1;
for (int i = 0; i < n; ++i) {
  int k = i;
  for (int j = i + 1; j < n; ++j)
    if (abs(a[j][i]) > abs(a[k][i])) k = j;
  if (abs(a[k][i]) < EPS) {
    det = 0;
    break;
  }
  swap(a[i], a[k]);
  if (i != k) det = -det;
  det *= a[i][i];
  for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
  for (int j = 0; j < n; ++j)
    if (j != i && abs(a[j][i]) > EPS)
      for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
}

cout << det;

• 无向图的生成树计数：矩阵树定理 在操作中一般直接计算上面这个行列式的值。