#include <iostream.h>
#include <math.h>

bool Luf_cp(double *a,int *cpp,int n,double eps=1.0e-8);      //作LU分解
void solve_lu_cpp(double *b,double *lu,int *cpp,int n,double eps=1.0e-8); //求AX=b的解

int main(int argc, char* argv[])
{
	//變量的數目
	int n; 
	//系數矩陣
	double *a; 
	//當前行主元變量序號
	int *cpp;
	//常數列
	double *b;  
	int i;

	cout<<"請輸入變量值"<<endl;
	cin>>n;
	if(n<1)  //非法變量數
	{
		cout<<"請輸入不小于1的數";
		return 0;
	}
	
	//分配空間
	a=new double[n*n];
	b=new double[n];
	cpp=new int[n];
	
	//讀入系數矩陣
	cout<<"請輸入系數矩陣"<<endl;
	for(i=0;i<n*n;i++)
		cin>>a[i];

	//讀入常數列
	cout<<"請輸入常量列"<<endl;
	for(i=0;i<n;i++)
		cin>>b[i];
	
	//初始化變量位置向量
	for(i=0;i<n;i++)
		cpp[i]=i;
	
	/*如果LU分解成功
	 *求解Ax=b*/
	if(Luf_cp(a,cpp,n))   
	{
		solve_lu_cpp(b,a,cpp,n); 

		cout<<"解答如下:"<<endl;
		for(i=0;i<n;i++)
			cout<<b[i]<<' ';
		cout<<endl;
	}
	else
		cout<<"sorry，不能用高斯方法解決!"<<endl;

	//輸出LU矩陣
	cout<<endl<<"LU矩陣如下:"<<endl;
	for(i=0;i<n*n;i++)
	{
		if(i%n==0)
			cout<<endl;
		cout<<a[i]<<'\t';
	}
	
	cout<<endl;
	//釋放存儲空間
	delete []a;
	delete []b;
	delete []cpp;

	return 0;
}

/*LU分解，返回職位是否成功，
 *若矩陣非奇異返回true，否則false
 *a輸入為系數矩陣，返回時為LU矩陣（左下角為響應的變換向量）
 *cpp為變換后各變量位置
 *n為變量數，eps為精度要求
 *絕對值小于eps的數近似為零*/
bool Luf_cp(double *a,int *cpp,int n,double eps)      
{
	int row,colume;

	//初始化變量位置向量
	for(row=0;row<n;row++)
		cpp[row]=row;

	for(colume=0;colume<n-1;colume++)
	{
		int maxrow=colume;

		//選擇當前列的最大值
		for(row=colume+1;row<n;row++)
			if(fabs(a[row*n+colume])>fabs(a[maxrow*n+colume]))
				maxrow=row;
		
		//如果整列為0，返回失敗
		if(a[maxrow*n+colume]==0)
			return false;

		//如果需要，交換最大值所在行與當前行
		if(maxrow!=colume)
		{
			double temp;
			int temp2;

			temp2=cpp[maxrow];
			cpp[maxrow]=cpp[colume];
			cpp[colume]=temp2;

			for(int col=0;col<n;col++) //交換
			{
				
				temp=a[maxrow*n+col];
				a[maxrow*n+col]=a[colume*n+col];
				a[colume*n+col]=temp;
			}
		}

		//變換一下各行
		for(row=colume+1;row<n;row++)
		{
			double percent;
			
			percent=a[row*n+colume]/a[colume*n+colume];
			a[row*n+colume]=percent;
			for(int col=colume+1;col<n;col++)
			{
				a[row*n+col]-=percent*a[colume*n+col];
				if(fabs(a[row*n+col])<eps)
					a[row*n+col]=0.0;
			}
		}
	}
	return true;
}

/*求解Ax=b
 *a輸入為系數矩陣，返回時為LU矩陣（左下角為響應的變換向量）
 *cpp為變換后各變量位置
 *n為變量數，eps為精度要求
 *絕對值小于eps的數近似為零*/
void solve_lu_cpp(double *b,double *lu,int *cpp,int n,double eps) 
{
	int row,colume;
	double *tb=new double[n];

	//重排常數列,tb為重排后數組
	for(row=0;row<n;row++)
		tb[row]=b[cpp[row]];

	//計算常數列
	for(row=1;row<n;row++)
	for(colume=0;colume<row;colume++)
		tb[row]-=lu[row*n+colume]*tb[colume];

	//從第n行到第1行回溯求解
	for(row=n-1;row>=0;row--)
	{
		for(colume=row+1;colume<n;colume++)
			tb[row]-=lu[row*n+colume]*tb[colume];

		tb[row]/=lu[row*n+row];
		if(fabs(tb[row])<eps)
			tb[row]=0;
	}

	//按x1--xn順序重排解向量,b中是結果
	for(row=0;row<n;row++)
		b[row]=tb[row];

	delete []tb;
}