#include <iostream.h>
#include <math.h>
#define MAX 10

void main()
{
	    int j;
		double zh;
	    double E[MAX][MAX];//定義單位矩陣
		for (int i=0;i<MAX;i++)
			for(j=0;j<MAX;j++)
			{
				if(i==j)E[i][j]=1;
				else
					E[i][j]=0;
			}
		double y[MAX],s1,s2,x[MAX][MAX];
		double a[MAX][MAX];
		int p[MAX]; // 記錄分解過程中的行交換!
		int k,N;
		cout<<"請輸入矩陣的階數:";
        cin>>N;
		cout<<"請輸入矩陣:"<<endl;
        for (i=0;i<N;i++)
			for (j=0;j<N;j++)
			{
			   cin>>a[i][j];
			}

		for(k=0;k<N;k++)
			p[k] = k+1;  // 初始化，相當于單位陣
//****** 開始LU分解 *************
		for(k=0;k<N-1;k++)
		{
			int s,t; //s存儲當前列絕對值最大的元素所在的行號,t為中間變量
			s=k;
			for(i=k+1;i<N;i++)
			{
				if(fabs(a[i][k])>fabs(a[s][k]))
					s=i;
			} // 選主元素
			if(s!=k)
			{
				double m;
				for(j=0;j<N;j++)// 注意這里要交換k、s兩行所有的元素
				{
					m = a[k][j];
					a[k][j] = a[s][j];
					a[s][j] = m;
					zh=E[k][j];E[k][j]=E[s][j];E[s][j]=zh;
				}  // 交換行
				t= p[s];
				p[s] = p[k];
				p[k] = t;  //記錄當前行交換
			} 

			if(a[k][k]==0)
			{
				cout << "存在為零的主元素！LU分解不能進行!" << endl;
				continue;
			}// 判定主元素為零，矩陣奇異，方程組解不唯一 

			for(i=k+1;i<N;i++)
			{
				a[i][k] = a[i][k] / a[k][k];  // 矩陣A的嚴格下三角部分存儲L矩陣的嚴格下三角部分 
				for(j=k+1;j<N;j++)
				{
					a[i][j] = a[i][j] - a[i][k] * a[k][j];
				}  //計算U矩陣     
			}
		}
		//下面是求解過程
		for(k=0;k<N;k++)
		{
            y[0]=E[k][0];
	     	for(i=1;i<N;i++)
			{   
		       s1=0;
		       for(j=0;j<i;j++)
                  s1=s1+a[i][j]*y[j];
		       y[i]=E[k][i]-s1;
			}   
	    	x[N-1][k]=y[N-1]/a[N-1][N-1];
    		for(i=N-2;i>=0;i--)
			{
		        s2=0;
		        for(j=i+1;j<N;j++)
			       s2+=a[i][j]*x[j][k];
	 	      x[i][k]=(y[i]-s2)/a[i][i];
			}	
		}
		for(i=0;i<N;i++)//輸出方陣
			for(j=0;j<N;j++)
			{
			   cout<<x[i][j]<<"   ";
			   if(j==N-1)cout<<endl;
			}
} //程序結束