#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 @
2009-08-20 19:27 Huicpc217 阅读(174) |
评论 (0) |
编辑 收藏