//	7.3  矩陣的LU分解

// 這個(gè)小程序應(yīng)該很容易看懂,所以不加詳細(xì)注釋了.
// 注釋見(jiàn)另一LU分解程序 

#include	<stdio.h>
#include    <stdlib.h>
#include	<conio.h>
#include	<math.h>

#define	MAX_n	100
#define PRECISION	0.0000001

void MatrixInput( float A[][MAX_n], int m, int n )
{
	int i, j;
	float ftmp;

	printf("\n===Begin input Matrix elements===\n");

	for( i = 1; i <= m; ++i )
	{
		printf( "Input_Line %d : ", i );
		for( j = 1; j <= n; ++j )
		{
			scanf("%f",&ftmp);
			A[i][j]=ftmp;
		}
	}
}

int LU_De_no_select( float A[][MAX_n], int n )
{
	int i,k,r;

	for( r = 1; r < n; ++r )
	{
		for( i = r + 1; i <= n; ++i )
		{
			if( fabs( A[r][r] ) < PRECISION ) return 1;
			
			for( k = 1; k < r; ++k )
			{
				A[i][r] -= A[i][k] * A[k][r];
			}

			A[i][r] /= A[r][r];
		}

		for( i = r + 1; i <= n; ++i )
		{
			for( k = 1; k <= r; ++k )
			{
				A[r+1][i] -= A[r+1][k] * A[k][i];
			}
		}
	}

	return 0;
}

void LowTriaMatrixOutput( float A[][MAX_n], int n )
{
	int i,j;

	for( i = 1; i <= n; ++i )
	{
		printf( "\n" );

		for( j = 1; j <= n; ++j )
		{
			if( j < i)
				printf( "%f\t", A[i][j] );
			else if( j == i )
				printf( "%f\t", 1.0 );
			else
				printf( "%f\t", 0.0 );
		}
	}
}

void UpTriaMatrixOutput( float A[][MAX_n], int n )
{
	int i,j;

	for( i = 1; i <= n; ++i )
	{
		printf( "\n" );

		for( j = 1; j <= n; ++j )
		{
			if( j >= i )
				printf( "%f\t", A[i][j] );
			else
				printf( "%f\t", 0.0 );
		}
	}
}

void main( )
{
	int n;
	float A[MAX_n][MAX_n];

	printf( "\nInput n=" );
	scanf( "%d", &n );
	if( n >= MAX_n - 1 )
	{
		printf( "\n\007n must <%d!", MAX_n );
		exit(0);
	}

	MatrixInput( A, n, n );

	if( LU_De_no_select( A, n ) )
	{
		printf( "\nLU Fault!" );
	}
	else
	{
		printf( "\nOutput L:" );
		LowTriaMatrixOutput( A, n );
		printf( "\nOutput U:" );
		UpTriaMatrixOutput( A, n );
	}

	printf( "\n\n\007Press any key to quit!\n" );
	getch();
}

/*
	運(yùn)行實(shí)例：

Input n=3

===Begin input Matrix elements===
Input_Line 1 : 2 -1 -1
Input_Line 2 : 3 4 -2
Input_Line 3 : 3 -2 4

Output L:
1.000000        0.000000        0.000000
1.500000        1.000000        0.000000
1.500000        -0.090909       1.000000
Output U:
2.000000        -1.000000       -1.000000
0.000000        5.500000        -0.500000
0.000000        0.000000        5.454545

Press any key to quit!

*/