//////////////////////////////////////////////////////////////////////
//	程序6.2   Gauss-Seidel 迭代

#include	<stdio.h>
#include	<stdlib.h>
#include	<conio.h>
#include	<math.h>

#define	MAX_n	100                //本程序能處理的方程最大階數
#define PRECISION	0.0000001      //迭代精度,相鄰二步解之差的最大范數是不是小于這一值，小于則迭代結束
#define MAX_Number	1000           //最大迭代次數,超過則算發(fā)散

//輸入n元向量到數組x中
void VectorInput(float x[],int n)
{
	int i;

	for(i=1;i<=n;++i)
	{
		printf("x[%d]=",i);
		scanf("%f",&x[i]);
	}
}

//輸入m X n 階陣 到二元數組 A[][] 中
void MatrixInput(float A[][MAX_n],int m,int n)
{
	int i,j;
	printf("\n===Begin input Matrix elements===\n");
	for(i=1;i<=m;++i)
	{
		printf("Input_Line %d : ",i);
		for(j=1;j<=n;++j)
			scanf("%f",&A[i][j]);
	}
}

//輸出n元向量數組x
void VectorOutput(float x[],int n)
{
	int i;
	for(i=1;i<=n;++i)
		printf("\nx[%d]=%f",i,x[i]);
}

//前后二步迭代之差的最大范數是不是小于精度要求
int IsSatisfyPricision(float x1[],float x2[],int n)
{
	int i;

	for(i=1;i<=n;++i)
		if(fabs(x1[i]-x2[i])>PRECISION) return 1;
	return 0;
}

//G-S迭代,n為階數，x為初始向量，A為方程組系數的擴展矩陣
//請與Jacobi迭代的程序比較，只有一點差別
int Jacobi_(float A[][MAX_n],float x[],int n)
{
	float x_former[MAX_n]; //放前一步結果向量
	int i,j,k;

	printf("\nInput vector x0:\n");
	VectorInput(x,n);     

	k=0;//當前迭代步
	do{
		//顯示上一步結果
		for(i=1;i<=n;++i)
		{
			printf("\nx[%d]=%f",i,x[i]);
			x_former[i]=x[i];
		}
		printf("\n");

		//迭代計算
		for(i=1;i<=n;++i)
		{
			x[i]=A[i][n+1];
			for(j=1;j<=n;++j)
				if(j!=i)x[i]-=A[i][j]*x[j];//Only defference:here x[j] replace x_former[j]
			if(fabs(A[i][i])>PRECISION)
				x[i]/=A[i][i];
			else
				return 1; //滿足精度，迭代結束
		}
		++k;
	}while(IsSatisfyPricision(x,x_former,n) && k<MAX_Number);

	if(k>=MAX_Number)
		return 1; //發(fā)散
	else
	{
		printf("\nG-S %d times!",k); //收斂，顯示迭代了多少步
		return 0;
	}
}

void main()
{
	int n;
	float A[MAX_n][MAX_n],x[MAX_n];

	//輸入階數
	printf("\nInput n=");
	scanf("%d",&n);
	if(n>=MAX_n-1) //不能超過最大充許階數，否則數組A放不下
	{
		printf("\n\007n must <%d!",MAX_n);
		exit(0);
	}

	//輸入方程組增廣（擴展）矩陣
	MatrixInput(A,n,n+1);

	//調用G-S迭代函數
	if(Jacobi_(A,x,n))
		printf("\nG-S Failed!");
	else
	{
		printf("\nOutput Solution:");
		VectorOutput(x,n); //輸出結果
	}
	printf("\n\n\007Press any key to quit!\n");
	getch();
}

/*
	運行實例：(注意：輸入的是方程組的增廣系數矩陣)

Input n=3

===Begin input Matrix elements===
Input_Line 1 : 2 -1 -1 4
Input_Line 2 : 3 4 -2 11
Input_Line 3 : 3 -2 4 11

Input vector x0:
x[1]=0
x[2]=0
x[3]=0
  .
  .
  .
x[1]=3.000000
x[2]=1.000000
x[3]=1.000000

G-S 21 times!
Output Solution:
x[1]=3.000000
x[2]=1.000000
x[3]=1.000000

Press any key to quit!
*/