高斯消元

来源:互联网 发布:618淘宝也买酒 编辑:程序博客网 时间:2024/04/30 03:15
bool equal(double source, double target){    if (abs(source - target) <=  0.000001)        return true;    else return false;}vector<double> Gauss(vector<vector<double> > A, vector<double> b,int &state){    //state -1 无解 0 唯一解 1 无穷解    //N个方程,M个未知数    int N = A.size(), M = A[0].size();    state = 0;    //处理上三角矩阵    for (size_t i = 0; i < M; i++)    {        bool find = false;        //令对角线上的元素 !=0        for (size_t j = i; j < N; j++)        {            if (A[j][i] != 0)            {                swap(A[i], A[j]);                find = true;                break;            }        }        //多解的情况        if (!find)        {            state = 1;            continue;        }        //将当前列第i行后的所有行都化为0(非上三角部分的化为0)        for (size_t j = i + 1; j < N; j++)        {            //double aji = A[j][i],aii= A[i][i];            double tmp = A[j][i] / A[i][i];            for (size_t k = 0; k < M; k++)                A[j][k] -= A[i][k] * tmp;                //A[j][k] -= A[i][k] * aji / aii;            b[j]-=b[i] * tmp;        }    }    //判断无解    for (size_t i = 0; i < N; i++)    {        bool allZero=true;        for (size_t j = 0; j < M; j++)        {            if (!equal(A[i][j],0))            {                allZero = false;                break;            }        }        if (allZero&&!equal(b[i],0))        {            state = -1;            return vector<double>();        }    }    //判断无穷解    if (state == 1) return vector<double>();    //求值过程    vector<double> result(M, 0);    for (int i = M - 1; i >= 0; i--)    {        for (size_t j = i + 1; j < M; j++)        {            b[i] -= A[i][j] * result[j];            A[i][j] = 0;        }        result[i] = b[i] / A[i][i];    }    return result;}
0 0