算法学习笔记1.1.2 高斯消元

任务

给一个n元一次方程组，求它们的解集。

说明

将方程组做成矩阵形式，再利用三种初等矩阵变换，得到上三角矩阵，最后回代得到解集。

接口

int solve(double a[][maxn], bool[], double ans[], const int& n);

• 复杂度：O(n^3)
• 输入：
  • a 方程组对应的矩阵
  • n 未知数个数
  • l,ans 存储解，l[]表示是否为自由元

代码

#include <iostream>
#include <cmath>
using namespace std;
const double EPS = 1e-9;
const int maxn = 1010;

inline int solve(double a[][maxn], bool l[], double ans[], const int& n) {
    int res = 0, r = 0;
    for (int i = 0; i < n; i ++)
        l[i] = false;
    for (int i = 0; i < n; i ++) {
        for (int j = r; j < n; j ++)
            if (fabs(a[j][i]) > EPS) {
                for (int k = i; k <= n; k ++)
                    swap(a[j][k], a[r][k]);
                break;
            }
        if (fabs(a[r][i]) < EPS) {
            res ++;
            continue;
        }
        for (int j = 0; j < n; j ++)
            if (j != r && fabs(a[j][i]) > EPS) {
                double tmp = a[j][i] / a[r][i];
                for (int k = i; k <= n; k ++)
                    a[j][k] -= tmp * a[r][k];
            }
        l[i] = true, r++;
    }
    for (int i = 0; i <n ; i ++)
        if (l[i])
            for (int j = 0; j < n; j ++)
                if (fabs(a[j][i]) > 0)
                    ans[i] = a[j][n] / a[j][i];
    return res;
}

使用范例

POJ1830