///////////////////////////////////////////////////////////////
//	§7  線性方程組的直接法
//
//	程序7.1  Gauss 消元法 — 不選主元

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

#define	MAX_n	100              //本程序能處理的方程最大階數
#define PRECISION	0.0000001    //主元是否小于這個數

//輸入m X n 階陣 到二元數組 A[][] 中
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;}
	}
}

//輸出 nXm 階矩陣A中的第k列,因為是輸出列，故m不用作為參數傳過來
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]);
}

//解n階上三解方程組,注意最后一列是方程組的右端向量
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;
}

//Gauss不選主元消去法，注意A的第n+1列是方程組的右端向量，1,...,n列是系數矩陣
int GaussElimination(float A[][MAX_n],int n)
{
	int i,j,k;
	
	//消元，最后成了上三角方程組
	for(i=1;i<n;++i)
	{
		if(fabs(A[i][i])<PRECISION) return 1; //主元太小，消元失敗
		for(j=i+1;j<=n;++j)
			for(k=i+1;k<=n+1;++k)
				A[j][k]-=A[i][k]*A[j][i]/A[i][i];
	}

	//解下三角方程組
	UpTriangle(A,n);
	return 0;
}

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

	//輸入階數到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(GaussElimination(A,n))
		printf("\nGauss Failed!");       //失敗
	else
	{
		printf("\nOutput Solution:");
		MatrixOneColumnOutput(A,n,n+1); //輸出解向量
	}
	printf("\n\n\007Press any key to quit!\n");
	getch();
}

/*
	運行實例：（注意，輸入為方程組的增廣矩陣）

Input n=3

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

Output Sulution:
x[1]=0.250000
x[2]=1.500000
x[3]=0.250000

Press any key to quit!
*/