高斯消元

来源:互联网 发布:淘宝以下影片允许出售 编辑:程序博客网 时间:2024/06/10 08:45
#include<iostream>#include<cstdio>#include<vector>#include<algorithm>#include<cstring>using namespace std;const double eps=1e-8;const int MAXN=1e3;typedef double Matrix[MAXN][MAXN];Matrix a;//m行n+1列,从0开始//增广矩阵,即a[i][n]是方程右边的常数//结束后a[i][n]是第i个未知数的值int dcmp(double x)//浮点数比较{    if(fabs(x)<eps) return 0;    return x<0?-1:1;}void Gauss(Matrix &a,int m,int n){    //变成行阶梯形矩阵    for(int i=0;i<m;++i)    {        //列选主元        int r=i;        for(int j=i+1;j<m;++j)        {            if(fabs(a[j][i])>fabs(a[r][i])) r=j;        }        if(r!=i) for(int j=0;j<=n;++j) swap(a[r][j],a[i][j]);        //与第i+1~m-1行进行消元        for(int k=i+1;k<m;++k)        {            double f=a[k][i]/a[i][i];            for(int j=i;j<=n;++j) a[k][j]-=f*a[i][j];        }    }}//n元线性方程组int solve(Matrix &a,int m,int n){    Gauss(a,m,n);    int r1=0;//系数矩阵的秩    int r2=0;//增广矩阵的秩    for(int i=0;i<m;++i)    {        bool zero=true;        for(int j=0;j<=n;j++)        {            if(dcmp(a[i][j]))            {                if(j<n) r1++;                r2++;                zero=false;                break;            }        }        if(zero||r1<r2) break;    }    if(r1<r2) return -1;//无解    else if(r1==r2&&r1<n) return 0;//有无穷多解    else    {        //回代        for(int i=n-1;i>=0;--i)        {            for(int j=i+1;j<n;++j)            {                a[i][n]-=a[j][n]*a[i][j];            }            a[i][n]/=a[i][i];        }        return 1;//有唯一解    }}int main(){    int n,m;    while(scanf("%d%d",&n,&m)!=EOF)    {        for(int i=0;i<m;++i)        {            for(int j=0;j<=n;++j)            {                scanf("%lf",&a[i][j]);            }        }        int ans=solve(a,m,n);        if(ans<0) printf("No solutions\n");        else if(ans==0) printf("Many solutions\n");        else        {            for(int i=0;i<n;++i)            {                printf("%d\n",(int)(a[i][n]+0.5));            }        }    }    return 0;}
原创粉丝点击