/***********************************************************/
/*                    Jacobi迭代法                         */
/***********************************************************/
#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#define n 3
double G[n][n];
////////////打印函數/////////////
void print()
{
	printf("G is:\n");
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{	
			printf("%8.3f,",G[i][j]);
			if(j==n-1) printf("\n");
		}
}
void print1(char ar[2],double *arr)
{
	printf("%s is:",ar);
	for(int j=0;j<n;j++)				
	{	
		if(j==0) printf("\t[");
		printf("%14.9f",*(arr+j));
		if(j!=n-1) printf(",");
		else printf("\t]");
	}
	printf("\n");
}
///////求數組的最大值////////////
double arraymax(double *array)
{
	double max=0;
	int maxnum=0;
	for(int kk=0;kk<n;kk++)
		if(max<fabs(*(array+kk))) {max=fabs(*(array+kk));maxnum=kk;}
	//printf("maxnum is %d,max is %8.6f\n",maxnum,max);
	return max;
}

void Jacobi(double A[n][n],double b[n],double ee)
{
	int i,j;
	int count=0;
	double temp=0;
	double Mid[n];
	double x[n]={0};
	double xx[n];
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if (A[i][i]==0) {printf("A[%d][%d]=0,error!",i,i);break;}
			if(i==j) G[i][j]=0;
			else G[i][j]=-A[i][j]/A[i][i];
		}
	print();//print array G[n][n];
	while(1)
	{
		for(i=0;i<n;i++)
		{	
			temp=0;
			for(j=0;j<n;j++)
				{	
				temp=temp+G[i][j]*x[j];
				}
			xx[i]=temp+b[i]/A[i][i];
			Mid[i]=xx[i]-x[i];
			//print1("x",x);
			//print1("xx",xx);
			//print1("Mid",Mid);
		}
		if(arraymax(Mid)<ee) break;
		else {
			for(i=0;i<n;i++)
				x[i]=xx[i];	
			count++;
		}
	}
	printf("Jocobi diedai count is %d\n",count);
	print1("the final x",x);
}
int main(int argc, char* argv[])
{
	double A1[n][n]={10,-2,-2,-2,10,-1,-1,-2,3};
    double b1[n]={1,0.5,1};
	double A[n][n]={5,2,1,-1,4,2,2,-5,10};
	double b[n]={-12,10,1};
	Jacobi(A,b,0.0001);
	Jacobi(A1,b1,0.00001);
	return 0;
}