/////////////////////////////////////
//	程序7.4  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;
		}
	}
}

// 不選主元LU分解
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;
}

// 下三角方程組求解
int LowTriangle_1( float L[][MAX_n], int n )
{
	int i,j;

	for( i = 1; i <= n; ++i )
//	{
//		if(fabs(L[i][i])<PRECISION)return 1;
		for( j = 1; j < i; ++j )
		{
			L[i][n+1] -= L[i][j] * L[j][n+1];
		}
//		L[i][n+1]/=L[i][i];
//	}

	return 0;
}

// 上三角方程組求解
int UpTriangle(float U[][MAX_n],int n)
{
	int i,j;

	for( i = n; i > 0; --i )
	{
		if( fabs( U[i][i] ) < PRECISION ) return 1;

		for( j = i + 1; j <= n; ++j )
		{
			U[i][n+1] -= U[i][j] * U[j][n+1];
		}

		U[i][n+1] /= U[i][i];
	}

	return 0;
}

// 矩陣第k列輸出, 用來輸出解向量
void MatrixOneColumnOutput( float A[][MAX_n], int n, int k )
{
	int i;

	for( i = 1; i <= n; ++i )
	{
		printf( "\nx[%d]=%f", i, A[i][k] );
	}
}

void main()
{
	int n;
	static 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 + 1 );

	if( LU_De_no_select( A, n ) )
	{
		printf("\nLU Fault!");
	}
	else
	{
		LowTriangle_1( A, n );
		
		if( UpTriangle( A, n ) )
		{
			printf( "\nLU Fault!" );
		}
		else
		{
			MatrixOneColumnOutput( A, n, n + 1 );
		}
//		printf("\nOutput L:");
//		LowTriaMatrixOutput(A,n);
//		printf("\nOutput U:");
//		UpTriaMatrixOutput(A,n);

	}

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

/*
	運行實例：(注意：輸入的是方程組的增廣系數矩陣)

Input n=3

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

x[1]=3.000000
x[2]=1.000000
x[3]=1.000000
*/