#include "stdio.h"
#include "stdlib.h"
#include "math.h"

int FunctionGauss(double a[][3],double b[],int n);//高斯消元法
int FunttionJacobi(double a[][3],double b[],double x[],double dbLimit,int nNum,int n);//雅克比迭代法



void main()
{
	//系數矩陣
	double dbM[3][3]=
	{{10,-1,-2},
	{-1,10,-2},
	{-1,-1,5}};
	
	
	double dbB[3]={7.2,8.3,4.2};

	//輸出解的結果
	if (FunctionGauss(dbM,dbB,3)==0)
	{
		printf("高斯消元法求解結果:\n");
		for(int i=1;i<=3;i++)
			printf("x(%d)=%e\n",i,dbB[i-1]);
	}
	else
	{
		printf("不能消元\n");
	}

	//雅克比迭代法求解

	double dbM1[3][3]=
	{{10,-1,-2},
	{-1,10,-2},
	{-1,-1,5}};
	
	
	double dbBj[3]={7.2,8.3,4.2};
	double dbIx[3]={0.0,0.0,0.0};
	printf("雅克比迭代法求解結果:\n");
	int nret=FunttionJacobi(dbM1,dbBj,dbIx,0.001,100,3);
	if (nret==1)
	{
		printf("雅克比迭代失效");
	}
	else if(nret==2)
	{
		printf("迭代不收斂");
	}

	getchar();
}

//函數名稱：FunctionGauss
//參數說明：a：系數矩陣，b：常數項，n：系數矩陣的維數
//返回值：1：消元失敗，0：消元成功
int FunctionGauss(double a[][3],double b[],int n)
{
	int k=0;
	int i=0;
	int j=0;
	//進行消元
	for(k=0;k<n-1;k++)
	{
	
		if (a[k][k]==0)
		{
			return 1;
		}
		for(i=k+1;i<=n-1;i++)
		{
			a[i][k]=a[i][k]/a[k][k]*-1.0;
			
			for(j=k+1;j<=n-1;j++)
			{
				a[i][j]+=a[i][k]*a[k][j];
			}
			b[i]+=a[i][k]*b[k];
		
		}
		
	}

	//回代
	b[n-1]/=a[n-1][n-1];

	double dbSum=0.0;
	for(k=n-2;k>=0;k--)
	{
		dbSum=0.0;
		for(j=k+1;j<=n-1;j++)
			dbSum+=a[k][j]*b[j];

		b[k]=(b[k]-dbSum)/a[k][k];
	}

	return 0;
}
//函數名稱：FunttionJacobi
//參數：a：系數矩陣，b：常數項，x0：初始值，dblimit：終止條件，nNum：最大迭代次數,n：矩陣維數
//返回值：0：求解完成，1：迭代失效，2：迭代不收斂
//
int FunttionJacobi(double a[][3],double b[],double x0[],double dbLimit,int nNum,int n)
{
	int k=0;
	int i=0;
	int j=0;
	double x[3]={0.0,0.0,0.0};
	double dbSum=0.0;
	double dbFabs=0.0;
	double dbMax=0.0;

	for(k=1;k<=nNum;k++)
	{
		for(i=0;i<n;i++)
		{

			if (a[i][i]==0)
				return 1;

			dbSum=0.0;
			for(j=0;j<n;j++)
			{
				if (j!=i)
				{
					dbSum+=a[i][j]*x0[j];
				}
			}

			x[i]=(b[i]-dbSum)/a[i][i];
		
		}
		dbMax=0.0;
		printf("第%d次迭代結果：\n",k);
		for(j=0;j<n;j++)
		{
			dbFabs=fabs(x[j]-x0[j]);
			x0[j]=x[j];
			printf("x(%d)=%e\n",j,x0[j]);
			if (dbFabs>dbMax)
				dbMax=dbFabs;
		}

		if (dbMax<dbLimit)
		{	printf("雅克比迭代結束");
			return 0;
		}
	
	}



	return 2;
}