浮点数高斯约当消元法模板

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define MAXN 400
#define EPS 1e-8
int n,m,row;
bool f[MAXN+10];
double a[MAXN+10][MAXN+10],x[MAXN+10];
template<class T>
void Read(T &x){
    char c;
    bool f=0;
    while(c=getchar(),c!=EOF){
        if(c=='-')
            f=1;
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            if(f)
                x=-x;
            return;
        }
    }
}
void read(){
    Read(n),Read(m);
    int i,j;
    for(i=1;i<=n;i++){
        for(j=1;j<=m;j++)
            Read(a[i][j]);
        Read(a[i][m+1]);
    }
}
void gauss_jordan(){
    int col,i,mxr,j;
    for(row=col=1;row<=n&&col<=m;row++,col++){
        mxr=row;
        for(i=row+1;i<=n;i++)
            if(fabs(a[i][col])>fabs(a[mxr][col]))
                mxr=i;
        if(mxr!=row)
            swap(a[row],a[mxr]);
        if(fabs(a[row][col])<EPS){
            row--;
            continue;
        }
        for(i=1;i<=n;i++)
            if(i!=row&&fabs(a[i][col])>EPS)
                for(j=m+1;j>=col;j--)
                    a[i][j]-=a[row][j]/a[row][col]*a[i][col];
    }
    row--;
}
void print(){
    int i,cnt,ff,j;
    if(row==n&&m==n){
        for(i=1;i<=n;i++)
            printf("x%d = %.4lf\n",i,a[i][m+1]/a[i][i]+EPS);
        return;
    }
    else{
        for(i=row+1;i<=n;i++)
            if(fabs(a[i][m+1])>EPS){
                puts("No solution");
                return;
            }
        printf("Multiple solution, free_num:%d\n",m-row);
        for(i=1;i<=row;i++){
            cnt=0;
            for(j=1;j<=m;j++)
                if(fabs(a[i][j])>EPS&&!f[j])
                    cnt++,ff=j;
            if(cnt>1)
                continue;
            x[ff]=a[i][m+1]/a[i][ff];
            f[ff]=1;
        }
        for(i=1;i<=m;i++)
            if(f[i])
                printf("x%d = %.4lf\n",i,x[i]+EPS);
            else
                printf("x%d not determined\n",i);
    }
}
int main()
{
    read();
    gauss_jordan();
    print();
}
posted @ 2016-01-27 23:13  outer_form  阅读(154)  评论(0编辑  收藏  举报