矩阵算法 高斯消元 行列式 矩阵的秩

今天学习一下矩阵的基本算法

高斯消元是解线性方程组的有力工具。

基本思想是通过将增广矩阵经过行初等变化变成简化阶梯形矩阵。

下面采用的是列主元高斯消元法，复杂度为O(n^3)。

很容易根据高斯消元法的过程得出行列式和秩的算法。

代码：

/*********************************************************
*            ------------------                          *
*   author AbyssalFish                                   *
**********************************************************/
#include<cstdio>
#include<iostream>
#include<string>
#include<cstring>
#include<queue>
#include<vector>
#include<stack>
#include<vector>
#include<map>
#include<set>
#include<algorithm>
#include<cmath>
//#include<bits/stdc++.h>
using namespace std;

typedef long long ll;


const double EPS = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;

//O(n^3)
vec gauss_jordan(const mat& A, const vec& b)
{
    int n = A.size();
    mat B(n,vec(n+1)); //Augment Matrix
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++) B[i][j] = A[i][j];
    for(int i = 0; i < n; i++) B[i][n] = b[i];

    for(int i = 0; i < n; i++){
        int piv = i; //取最大以便判断无解或者无穷多解
        for(int j = i; j < n; j++){
            if(abs(B[j][i]) > abs(B[piv][i])) piv = j;
        }
        if(i != piv) swap(B[i],B[piv]);
        if(abs(B[i][i]) < EPS) return vec();

        //假想把系数变成1，只计算对后面有影响的部分
        for(int j = n; j > i; j--) B[i][j] /= B[i][i];
        for(int j = 0; j < n; j++) if(i != j) {
            for(int k = i+1; k <= n; k++) B[j][k] -= B[j][i]*B[i][k];
        }
    }
    vec x(n);
    for(int i = 0; i < n; i++) x[i] = B[i][n];
    return x;
}


double determinant(const mat& A)
{
    int n = A.size();
    mat B = A;
    double det = 1;
    int sign = 0;
    for(int i = 0; i < n; i++){
        int piv = i;
        for(int j = i; j < n; j++){
            if(abs(B[j][i]) > abs(B[piv][i])) piv = j;
        }
        if(i != piv) swap(B[i],B[piv]), sign ^= 1;
        if(abs(B[i][i]) < EPS) return 0;
        det *= B[i][i];
        for(int j = i+1; j < n; j++) B[i][j] /= B[i][i];
        for(int j = i+1; j < n; j++) {
            for(int k = i+1; k < n; k++) B[j][k] -= B[j][i]*B[i][k];
        }
    }
    return sign?-det:det;
}

int rank_of_mat(const mat& A)
{
    int n = A.size();
    mat B = A;
    for(int i = 0; i < n; i++){
        int piv = i;
        for(int j = i; j < n; j++){
            if(abs(B[j][i]) > abs(B[piv][i])) piv = j;
        }
        if(i != piv) swap(B[i],B[piv]);
        if(abs(B[i][i]) < EPS) return i;
        for(int j = n; --j > i;) B[i][j] /= B[i][i];
        for(int j = i+1; j < n; j++) {
            for(int k = i+1; k < n; k++) B[j][k] -= B[j][i]*B[i][k];
        }
    }
    return n;
}

double read(){ double t; scanf("%lf",&t); return t; }

//#define LOCAL
int main()
{
#ifdef LOCAL
    freopen("in.txt","r",stdin);
#endif
    vec alpha;
    mat A;
    int n;
    puts("input a matrix");
    while(~scanf("%d",&n) && n <= 0) puts("input invalid");
    A.resize(n);
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            A[i].push_back(read());
        }
    }
    puts("input a vector");
    for(int i = 0; i < n; i++) alpha.push_back(read());

    puts("
solution");
    vec sol = gauss_jordan(A,alpha);
    if(sol.size()){
        for(int i = 0; i < n; i++)
            printf("%lf%c", sol[i], i!=n-1?' ':'
');
    }else {
        puts("not unique");
    }

    puts("
determinant");
    printf("%lf
", determinant(A));

    puts("
rank");
    printf("%d
", rank_of_mat(A));


    return 0;
}

测试样例

实际使用中如果同时需要这些信息只计算一遍B矩阵，从而提高效率。