高斯消元

Huicpc217 — Thu, 20 Aug 2009 11:27:00 GMT

#include<iostream>
#include<math.h>
using namespace std;

const double eps = 1e-10;

inline int sig( double x){
    return (x > eps) - (x < - eps);
}

inline void swap( double &a,double &b){
    double t = a; a = b ; b = t;
}

int lup_decomposition( double **a , int n , int pi[])
{
    // n = rows[ A ];
    int i, j, k, k1;
    double p;
    for ( i = 0 ; i < n ; i ++ )
        pi[i] = i;// 置换

    for ( k = 0 ; k < n ; k ++ ){
        p = 0;
        for ( i = k ; i < n ; i ++ )
        {
            if( fabs( a[i][k] ) > p )
            {
                p = fabs( a[i][k] );
                k1 = i;
            }
        }
        if( sig( p ) == 0 ){
            return 0 ;// error
        }
        swap( pi[ k ] , pi[ k1 ] );
        for ( i = 0 ; i < n ; i ++ )
            swap( a[k][i], a[k1][i] );
        for ( i = k + 1 ; i < n ; i ++ )
        {
            a[i][k] = a[i][k] / a[k][k];
            for ( j = k + 1 ; j < n ; j ++ )
                a[i][j] = a[i][j] - a[i][k] * a[k][j] ;
        }
    }
    return 1;
}

void lup_solve(double * x,double * y, double **L,double **U,int pi[], double *b, int n)
{
    // n = rows[ L ]
    int i, j;
    double sum;
    for ( i = 0 ; i < n ; i ++ ){
        sum = 0;
        for ( j = 0 ; j < i ; j ++ )
            sum += L[i][j] * y[ j ];
        y[i] = b[ pi[i] ] - sum;
    }
    for ( i = n - 1 ; i >= 0 ; i -- ){
        sum = 0;
        for ( j = i + 1 ; j < n ; j ++ )
            sum += U[i][j] * x[ j ];
        x[i] = (y[i] - sum)/U[i][i];
    }
}

int main()
{
    int i, j, k;
    int n;
    int z = 0;
    while( scanf("%d",&n) != EOF )
    {
        if( z ) printf("\n");
        z = 1;
        int * pi = new int[ n ];
        double * b = new double [ n ], * x = new double [n] , * y = new double [n];
        double ** a = new double *[n];
        for ( i = 0 ; i < n ; i ++ )
            a[i] = new double [n];

        for ( i = 0 ; i < n ; i ++ ){
        for ( j = 0 ; j < n ; j ++ )
            scanf("%lf",&a[i][j]);
            scanf("%lf",&b[i]);
        }
        //for ( i = 0 ; i < n ; i ++ )
          //  scanf("%lf",&b[i]);

        int  flag ;
        flag = lup_decomposition( a, n, pi );
        if( flag )
        {
            lup_solve( x, y, a, a, pi, b, n );
            for ( i = 0 ; i < n ; i ++ )
                cout<<x[i]<<' ';
            cout<<endl;
        }
        else
            cout<<"No solution."<<endl;

    }
    return 0;
}
/*

3
1 2 0
3 4 4
5 6 3

3 7 8

*/

Huicpc217 2009-08-20 19:27 发表评论

test

Huicpc217 — Wed, 19 Aug 2009 07:18:00 GMT

justs a test

C++博客-^^想你

高斯消元

test