//  cpp  t2_1b源代碼  上機練習題2_1b"
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
const int N=7;
static double aa[N][N+1]=
{
{3,-5,6,4,-2,-3,8,11},
{1,1,-9,15,1,-9,2,2},
{2,-1,7,5,-1,6,11,29},
{-1,1,3,2,7,-1,-2,9},
{4,3,1,-7,2,1,1,5},      
{2,9,-8,11,-1,-4,-1,8},
{7,2,-1,2,7,-1,9,25}
},aa2[N][N+1];
double x[N+1];
int main()
{
	int i,j,k,u;
	double m;
    void ShowResult(double aa[N][N+1],double x[N+1]) ; //函數聲明
	for(i=0;i<N;i++)             ///數據備份
		for(j=0;j<N+1;j++)
			aa2[i][j]=aa[i][j];
	cout<<"上機練習題2_1b"<<endl;
	cout<<"增廣矩陣(A|b)為:"<<endl;
	for(i=0;i<N;i++)                       ///輸出矩陣A|b
	{		
		for(j=0;j<N+1;j++)
			cout<<setw(9)<<aa[i][j];
		cout<<endl;
	}
	////Gauss消去過程
	for(k=0;k<N-1;k++)                     
	{
		if(fabs(aa[k][k])<1E-6)             ////主元不能為0
		{	cout<<"主元為0，計算不能進行！！！"<<endl;
			break;
		}
		for(i=k+1;i<N;i++)
		{
			aa[i][k]=aa[i][k]/aa[k][k];		
			for(j=k+1;j<=N;j++)
			{
				aa[i][j]=aa[i][j]-aa[i][k]*aa[k][j];
			}
		}
	}
	cout<<"Gauss消去后的增廣矩陣(A|b)為:"<<endl;
	ShowResult(aa,x); 

     //////列主元Gauss消去法
	for(k=0;k<N-1;k++)                     ////消去過程
	{
		u=k;
		for(i=k;i<N;i++)
			if(fabs(aa2[i][k])>fabs(aa2[u][k]))
				u=i;
		if(fabs(aa2[u][k])<1E-6)             ////主元不能為0
		{	cout<<"主元為0，計算不能進行！！！"<<endl;
			break;
		}  
		if(u!=k)
			for(j=0;j<N+1;j++)        ////交換k,r兩行
			{
				m=aa2[k][j];aa2[k][j]=aa2[u][j];aa2[u][j]=m;
			}
		
		for(i=k+1;i<N;i++)       ///消元計算
		{
			aa2[i][k]=aa2[i][k]/aa2[k][k];		
			for(j=k+1;j<=N;j++)
			{
				aa2[i][j]=aa2[i][j]-aa2[i][k]*aa2[k][j];
			}
		}	
	}
	cout<<"列主元Gauss消去后的增廣矩陣(A|b)為:"<<endl;
	ShowResult(aa2,x); 
	return 1;
}

void ShowResult(double aa[N][N+1],double x[N+1])   ////回代和顯示結果
{
	int i,j,k;
	double ss;
	for(i=0;i<N;i++)                       
	{		
		for(j=0;j<N+1;j++)
			cout<<setw(9)<<aa[i][j];
		cout<<endl;
	}
	                                       /////回代過程
	x[N-1]=aa[N-1][N]/aa[N-1][N-1];
	for(k=N-2;k>=0;k--)
	{
		ss=aa[k][N];
		for(j=k+1;j<N;j++)
			ss=ss-aa[k][j]*x[j];
		x[k]=ss/aa[k][k];
	}
	                                       ////輸出結果
	cout<<"計算結果為："<<endl;
	for(k=0;k<N;k++)
	{
		cout<<"x["<<k+1<<"]="<<x[k]<<"\t";
	}
	cout<<endl;
}