////////////////////////////////////////////////////////////////////////////////////////////////
#include"myhead.h"
////////////////////////////////////////////////////////////////////////////////////////////////
void exchange(int r,int k);
double max(int k);
double max(double a[],int n);
int flag,n;
double A[20][21],ark;
///////////////////////////////////////////////////////////////////////////////////////////////
void fun_2()
{
	 system("cls");//執行系統清屏命令
	 cout<<"	 ******************************************************************"<<endl
		 <<"	 *                        求解線性方程組的解                      *"<<endl
		 <<"	 ******************************************************************"<<endl
		 <<"	 *                       1、高斯列主元消去法：                    *"<<endl
	 	 <<"	 *                       2、矩陣的LU分解法：                      *"<<endl
		 <<"	 *                       3、雅可比迭代法：                        *"<<endl
		 <<"	 ******************************************************************"<<endl;
	 cout<<"	 請選擇一個數字(1-3):";
	 int ch=0;
	 cin>>ch;

	 if(ch==1)
		 g();
	 if(ch==2)
		 LU();
	 if(ch==3)
		 y();
	 system("pause");
}
//////////////////////////////////////////////////////////////////////////////////////////////
void g()
{
	double x[21];      //此數組用于存放方程解
	int k,i,j;
	printf("	 輸入方程組的維數:n=");
	scanf("%d",&n);
	printf("	 現在輸入系數矩陣A和向量b:");
	for(i=1;i<=n;i++)
	{
		printf("\n	 請輸入a%d1--a%d%d系數和向量b%d:",i,i,n,i);  //實現將每一行中的系數和向量一次性輸入， 
                                                               //數之間用空格格開，輸完后回車確定   
		for(j=1;j<=n+1;j++)     //將剛才輸入的數存入數組
		scanf("%f",&A[i][j]);
   }

   for(k=1;k<=n-1;k++)                       
   {
	   ark=max(k);
	   if(ark==0)                //判斷方程是否為線性方程，即是否合法
	   {
		   printf("\n	 此方程組不合法!");
		   return;
	   }
	   else if(flag!=k)
				exchange(flag,k);
		for(i=k+1;i<=n;i++)
			for(j=k+1;j<=n+1;j++)
				A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];
   }

   x[n]=A[n][n+1]/A[n][n];

   for( k=n-1;k>=1;k--)
   {
	   double me=0;
	   for(j=k+1;j<=n;j++)
	   {
		   me=me+A[k][j]*x[j];
		}
       x[k]=(A[k][n+1]-me)/A[k][k];
   }
   for(i=1;i<=n;i++)
   {
	   printf("\n x%d=%f",i,x[i]);
   }
   system("pause");
}

///////////////////////////////////////////////////////////////////////////////
void exchange(int r,int k)         //交換行的矩函數
{
	int i;
	for(i=1;i<=n+1;i++)
		A[0][i]=A[r][i];
	for(i=1;i<=n+1;i++)
		A[r][i]=A[k][i];
	for(i=1;i<=n+1;i++)
		A[k][i]=A[0][i];
}

////////////////////////////////////////////////////////////////////////////////
double max(int k)           //比校系數大小的函數
{
	int i;
	double  temp=0;
	for(i=k;i<=n;i++)
    if(fabs(A[i][k])>temp)
    {	
		temp=fabs(A[i][k]);
		flag=i;
    }
	return temp;
}
/////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////
void  y()                                 //雅可比函數
{
	double a[10][10],b[10],x0[10],x[10],y[10];
	int n;
	int i,j;
	double  sum,e;
	cout<<"	 請輸入方程組的階數和誤差限：";
	cin>>n>>e;
	cout<<"	 請輸入系數矩陣：";
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			cin>>a[i][j];
	cout<<"	 請輸入B矩陣：";
    for(i=0;i<n;i++)
		cin>>b[i];
	cout<<"	 請輸入X的初值："<<endl;
    for(i=0;i<n;i++)
		cin>>x0[i];
	while(1)
		{	for(i=0;i<n;i++)
				{  sum=0;
					for(j=0;j<n;j++)
						sum=a[i][j]*x0[j]+sum;
					x[i]=x0[i]+(b[i]-sum)/a[i][i];
					y[i]=fabs(x[i]-x0[i]);
				}
			if(max(y,n)<=e)break;
			else 
				for(i=0;i<n;i++)
					x0[i]=x[i];
		}
	cout<<"	 方程組的根為："<<endl;
    for(i=0;i<n;i++)
		cout<<x[i]<<endl;
}
///////////////////////////////////////////////////////////////////////////////////////////
double  max(double a[],int n)
{ 
	double m;
	m=a[0];
	for(int i=0;i<n;i++)
	  if(a[i]>m)
			 m=a[i];
	return m;
}
///////////////////////////////////////////////////////////////////////////////////////////////
//LU分解
int  LU()
{ 
	double a[10][10],b[10],l[10][10],u[10][10],x[10],y[10];
	int n,i,j,k;
	double a1,b1,c,d;
	cout<<"	 請輸入方程組的階數：";
	cin>>n;
	cout<<"	 請輸入系數矩陣"<<endl;
	for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
		{cin>>a[i][j];
		    l[i][j]=0;
			u[i][j]=0;}
    cout<<"	 請輸入B矩陣"<<endl;
    for(j=1;j<=n;j++)
		cin>>b[j];
    for(i=1;i<=n;i++)
	{u[1][i]=a[1][i];
	l[i][i]=1;}
    for(i=2;i<=n;i++)
		 l[i][1]=a[i][1]/u[1][1];
    for(i=2;i<=n;i++)
		{
		  for(j=i;j<=n;j++)
		  { b1=0;
			for(k=1;k<=i-1;k++)
				b1=b1+l[i][k]*u[k][j];
			u[i][j]=a[i][j]-b1;
		  }
		  for(j=i+1;j<=n;j++)
		  { a1=0;
	          for(k=1;k<=i-1;k++)
				a1=a1+l[j][k]*u[k][i];
			l[j][i]=(a[j][i]-a1)/u[i][i];
		  }
		}
	cout<<"	 產生的L矩陣為："<<endl;
     for(i=1;i<=n;i++)
		{for(j=1;j<=n;j++)
			cout<<l[i][j]<<"   ";
	     cout<<endl;
		}
	cout<<"	 產生的U矩陣為："<<endl;
     for(i=1;i<=n;i++)
		{for(j=1;j<=n;j++)
			cout<<u[i][j]<<"   ";
	     cout<<endl;
		}
	 y[1]=b[1];
	 for(i=2;i<=n;i++)
		{c=0;
		 for(k=1;k<=i-1;k++)
		  c=c+l[i][k]*y[k];
	      y[i]=b[i]-c;
		}
	 x[n]=y[n]/u[n][n];
	 for(i=n-1;i>=1;i--)
		{d=0;
	      for(j=i+1;j<=n;j++)
			  d=u[i][j]*x[j]+d;
		  x[i]=(y[i]-d)/u[i][i];
		}
	 cout<<"	 方程組的根為："<<endl;
	 for(i=1;i<=n;i++)
	 {cout<<'x'<<i<<'=';
	 cout<<x[i]<<endl;}
 return 0;
}