/*-矩陣的三角分解LU分解法-*/
#include<stdio.h>
#define N 10 /*-最多的未知數個數-*/
void main()
{
	int i,k,r,n,uk,lk;
	double a[N+1][N+1],b[N+1],y[N+1],x[N+1],l[N*(N+1)/2+1],u[N*(N+1)/2+1];
	double sum;
	printf("Input n which is lower 10:");
	scanf("%d",&n);
    
	/*-輸入系數矩陣A-*/
	printf("Input Matrix of A:\n");
	for(i=1;i<=n;i++)
		for(k=1;k<=n;k++)
			scanf("%lf",&a[i][k]);

	/*-輸入矩陣b-*/
	printf("Input Matrix of b:\n");
	for(i=1;i<=n;i++)
		scanf("%lf",&b[i]);
	
	/*-計算U的第一行元素-*/
	for(i=1;i<=n;i++)
	{
		uk=i*(i-1)/2+1;    /*-U中元素對應壓縮矩陣的位置-*/
		u[uk]=a[1][i];     /*-U的第一行元素-*/
	}

	/*-計算L的第一列元素-*/
	for(i=2;i<=n;i++)
	{
		lk=i*(i-1)/2+1;     /*-L中元素對應壓縮矩陣的位置-*/ 
		l[lk]=a[i][1]/u[1]; /*-L的第一行元素-*/ 
	}

	/*-對于r=2,3,4....,n進行計算-*/
	for(r=2;r<=n;r++)
	{
		/*-計算U的第r行元素-*/
		for(i=r;i<=n;i++)
		{
			sum=0;
			for(k=1;k<=r-1;k++)
			{
				uk=i*(i-1)/2+k;/*-轉成壓縮位置-*/
                lk=r*(r-1)/2+k;
				sum+=u[uk]*l[lk];
			}
			uk=i*(i-1)/2+r;
			u[uk]=a[r][i]-sum;
		}
		/*-計算L的第r列元素(r!=n)-*/
		if(r!=n)
		{
			for(i=r+1;i<=n;i++)
			{
				sum=0;
				for(k=1;k<=r-1;k++)
				{
					uk=r*(r-1)/2+k;
					lk=i*(i-1)/2+k;
					sum+=u[uk]*l[lk];
				}
				lk=i*(i-1)/2+r;
				uk=r*(r-1)/2+r;
				l[lk]=(a[i][r]-sum)/u[uk];
			}
		}/*-if()end-*/
	}/*-for(r)end-*/

	/*-求解Ly=b,Ux=y公式-*/
	y[1]=b[1];
	for(i=2;i<=n;i++)
	{
		sum=0;
		for(k=1;k<=i-1;k++)
		{
			lk=i*(i-1)/2+k;
			sum+=l[lk]*y[k];
		}
		y[i]=b[i]-sum;
	}
    x[n]=y[n]/u[n*(n-1)/2+n];
	for(i=n-1;i>=1;i--)
	{
		sum=0;
		for(k=i+1;k<=n;k++)
		{
			uk=k*(k-1)/2+i;
			sum+=u[uk]*x[k];
		}
		x[i]=(y[i]-sum)/u[i*(i-1)/2+i];
	}
	/*-輸出x1,x2,x3.....,xn-*/
	for(i=1;i<=n;i++)
		printf("x%d=%.8lf\n",i,x[i]);
}