高斯消元模板

最新推荐文章于 2022-04-12 19:24:58 发布

lshacm

最新推荐文章于 2022-04-12 19:24:58 发布

阅读量250

点赞数 1

分类专栏：矩阵与线性方程组

本文链接：https://blog.csdn.net/u013514182/article/details/52103540

版权

矩阵与线性方程组专栏收录该内容

7 篇文章 0 订阅

订阅专栏

浮点数高斯消元：

bool Free[maxm];

int Gauss(int n, int m, double A[][maxm]){//A[i][0..m-1]系数矩阵 A[i][m]:常数矩阵
    int row = 0;
    for(int col = 0; row < n && col < m; row ++, col ++){
        int mxr = row;
        for(int i = row + 1; i < n; i ++){
            if(fabs(A[i][col]) > fabs(A[mxr][col])) mxr = i;
        }
        if(mxr != row){
            for(int j = 0; j <= m; j ++) swap(A[row][j], A[mxr][j]);
        }

        if(fabs(A[row][col]) < eps) { row --; continue; }

        //将第col列出了第row行外全都消为0
        for(int i = 0; i < n; i ++){
            if(i == row) continue;
            double xi = A[i][col] / A[row][col];
            for(int j = m ;j >= col; j --){
                A[i][j] -= xi * A[row][j];
            }
        }
    }

    //判断解的情况

    //1.有唯一解
    if(row == n && m == n){
        for(int i = 0; i < n; i ++)
            x[i] = A[i][m] / A[i][i];
        return 0;
    }
    else{
        //2.无解

        for(int i = row; i < n; i ++)
            if(fabs(A[i][m]) > eps){
                return -1;
            }

        //3.多组解

        memset(Free, true, sizeof(Free));
        for(int i = row-1; i >= 0; i --){
            int cnt = 0, id = -1;
            //判断一行内自由变元的个数
            for(int j = 0; j < m; j ++){
                if(fabs(A[i][j]) > eps && Free[j]) { cnt++; id = j; } //该eps不可用0替代，否则会出问题
            }
            if(!cnt || cnt > 1) continue;

            //一行只有一个自由变元，可解
            double tmp = A[i][m];
            for(int j = 0; j < m; j ++){
                if(j != id && fabs(A[i][j]) > eps) tmp -= A[i][j] * x[j]; 
            }
            x[id] = tmp / A[i][id];
            Free[id] = 0;
        }
        return m - row;//返回自由变元的个数
    }
}