#include <stdio.h>
#include <math.h>
#define MAX_N 20

main()
{
	int n;
	int i,j,k;
	int mi;
	double mx,tmp;
	static double a[MAX_N][MAX_N],b[MAX_N],x[MAX_N];
    printf("實驗五  列主元高斯消元法解 Ax=b\n\n");
	printf("\n請輸入方程組的階數N：");
	scanf("%d",&n);
	if (n>MAX_N)
	{
		printf("你輸入的階數大于程序計算的范圍（20階），請重新輸入。\n");
	}
	else if (n<=0)
	{
		printf("階數N不能小于1，請重新輸入。（1到20）\n",MAX_N);
	}
	else
	{
	printf("請輸入方程組的系數（矩陣A）（系數間以空格鍵隔開）: \n");
	for (i=0;i<n;i++)
		for(j=0;j<n;j++)
		    scanf("%lf",&a[i][j]);
	printf("請輸入方程的右端項（矩陣b）%d: \n",n-1);
	for(i=0;i<n;i++) scanf("%lf",&b[i]);
	for(i=0;i<n-1;i++)
	{
		for(j=i+1,mi=i,mx=fabs(a[i][i]);j<n;j++)
		if (fabs(a[j][i])>mx)
		{
			mi=j;
			mx=fabs(a[j][i]);
		}
		if (i<mi)
		{
			tmp=b[i];b[i]=b[mi];b[mi]=tmp;
			for(j=i;j<n;j++)
			{
				tmp=a[i][j];
				a[i][j]=a[mi][j];
				a[mi][j]=tmp;
			}
		}
		for(j=i+1;j<n;j++)
		{
			tmp=-a[j][i]/a[i][i];
			b[j]=b[j]+b[i]*tmp;
			for(k=i;k<n;k++) a[j][k]=a[j][k]+a[i][k]*tmp;
		}
	}
	x[n-1]=b[n-1]/a[n-1][n-1];
	for(i=n-2;i>=0;i--)
	{
		x[i]=b[i];
		for(j=i+1;j<n;j++)
		x[i]=x[i]-a[i][j]*x[j];
	    x[i]=x[i]/a[i][i];		
	}
	printf("計算結果：\n");
	for (i=0;i<n;i++) printf("x[%d]=%f\n",i,x[i]);
	}
	getchar();
	getchar();
}