posts - 2,comments - 1,trackbacks - 0
#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

*/

posted on 2009-08-20 19:27 Huicpc217 阅读(168) 评论(0)  编辑 收藏 引用

只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理