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

int LU_De_no_select(float A[][MAX_n],int n)
{
	int i,j,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;
}

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;
	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 Failt!");
	else
	{
		LowTriangle_1(A,n);
		if(UpTriangle(A,n))
			printf("\nLU Failt!");
		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
*/