雕刻时光

just do it……nothing impossible
  博客园  :: 首页  :: 新随笔  :: 联系 :: 订阅 订阅  :: 管理

高斯消元求解线性方程组

Posted on 2013-10-23 19:27  huhuuu  阅读(739)  评论(2编辑  收藏  举报

http://202.120.80.191/problem.php?problemid=1040

裸题 

Description 

应用高斯消元法求解n*n的线性方程组Ax=b,其中A为系数矩阵。
数据保证有唯一解。

Input 

第1行为一个整数n(0<n<=20),表示是n*n的的线性方程组。
接下去的n行表示了系数矩阵A,每行有n个整数。
再接下去的n行表示了b,每行只有一个整数。

Output 

输出有n行,每行有1个小数(精确到0.01),表示方程组的解。

Sample Input 

3
1 2 3
2 4 5
3 1 2
4
3
1

Sample Output 

-1.40
-4.80
5.00

ps:手写了个高斯消元 AC,超爽 ,注意要对double 0的情况特殊处理

#include<stdio.h>
#include<math.h>
#include<iostream>
using namespace std;
#define eps 1e-8  
#define zero(a) fabs(a)<eps

double map[29][29];
double ret[29];

void mat(int n){
    int i,j,k,ti;
    for(k=1;k<n;k++){
        int firsti=-1;
        double map_i_k;
        for(i=k;i<=n;i++){
            if(firsti==-1){
                if(zero(map[i][k])){
                    for(ti=i+1;ti<=n;ti++){
                        if(!zero(map[ti][k]))break;
                    }
                    if(ti==n+1)break;
                    for(j=1;j<=(n+1);j++){
                        swap(map[i][j],map[ti][j]);
                    }
                }
                firsti=i;
            }

            if(zero(map[i][k]))continue;
            map_i_k=map[i][k];
            for(j=1;j<=(n+1);j++){
                map[i][j]/=map_i_k;
            }
            if(firsti==-1)continue;
            if(firsti==i)continue;
            for(j=1;j<=(n+1);j++){
                map[i][j]-=map[firsti][j];
            }
        }
    }

    //test
    //for(i=1;i<=n;i++)ret[i]=0;
    //for(i=1;i<=n;i++){
    //    for(j=1;j<=n;j++){
    //        if(!zero(map[i][j]))break;
    //    }
    //    if(j==n+1)continue;
    //    ret[j]=map[i][n+1]/map[i][j];
    //}
    //for(j=1;j<=n;j++){
    //    printf("%.2lf\n",ret[j]);
    //}

    //for(i=1;i<=n;i++){
    //    for(j=1;j<=(n+1);j++){
    //        printf("%.8lf ",map[i][j]);
    //    }printf("\n");
    //}
    //
    for(i=n;i>=1;i--){
        for(j=n;j>i;j--){
            map[i][n+1]-=map[i][j]*ret[j];
        }
        ret[i]=map[i][n+1]*1.0/map[i][i];
    }
}

int main(){
    int n;
    while(scanf("%d",&n)!=EOF){
        int i,j;
        for(i=1;i<=n;i++){
            for(j=1;j<=n;j++){
                scanf("%lf",&map[i][j]);
            }
        }
        for(i=1;i<=n;i++){
            scanf("%lf",&map[i][n+1]);
        }
        mat(n);
        for(i=1;i<=n;i++){
            if(zero(ret[i]))printf("0.00\n");
            else
            printf("%.2lf\n",ret[i]);
        }
    }

    return 0;
}
View Code