#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
//------------------------------
//Data[i][j];
//i - Stroka  (Row)
//j - Stolbec (Col)
//------------------------------
//------------------------------
//Boolean Type
#define MINVAL 0.0001f
typedef enum
{
	FALSE=0,
	TRUE
} BOOL;
//------------------------------
//Matrix Type
struct QMATRIX
{
	float **Data;
	int n;
};
//------------------------------
//Vector Type
typedef float *QVECTOR;
//------------------------------
#define  AllocQVector(n) ((QVECTOR)malloc(sizeof(float)*(n))) //Allocate Vector
#define  CopyQVector(v1,v2,n) (memcpy(v1,v2,sizeof(float)*(n))) //Copy From v2 To v1
#define  FreeQVector(Vector) free(Vector) //Delete QVECTOR
//------------------------------
//Prototypes
float DeltaQVectors(QVECTOR v1,QVECTOR v2,int n);
//------------------------------
void FreeQMatrix(QMATRIX *m);
BOOL LoadQMatrix(char *FileName,QMATRIX *m);
BOOL PrintQMatrix(QMATRIX *m);
BOOL QSwapRows(QMATRIX *m,int i,int j);
BOOL QIterateRows(QMATRIX *m,int From);
BOOL TransformQMatrix(QMATRIX *m);
void QGetNewX(QMATRIX *m,QVECTOR NewX);
QVECTOR QGaussZ(QMATRIX *m,float e);
BOOL QGaussDown(QMATRIX *m,int from);
BOOL QGaussUp(QMATRIX *m);
BOOL QGauss(QMATRIX *m);

void PrintQResult(QMATRIX *m,QVECTOR Result);
//------------------------------
//Free Matrix


void FreeQMatrix(QMATRIX *m)
{
	int i;

	if(m==NULL)
		return;
	m->n=0;
	if(m->Data==NULL)
		return;
	for(i=0;i<m->n;i++)
		if(m->Data[i]!=NULL)
			free(m->Data[i]);
	free(m->Data);
}
//------------------------------
//Load Matrix From File
BOOL LoadQMatrix(char *FileName,QMATRIX *m)
{
	FILE *f;
	int i,j;

	//Free Last Matrix
	FreeQMatrix(m);

	//Open File
	f=fopen(FileName,"rt");
	if(f==NULL)
		return FALSE;

	//Scan Size
	if(fscanf(f,"%i",&(m->n))!=1)
	{
		m->n=0;
		fclose(f);
		return FALSE;
	}

	//Allocate Memory
	m->Data=(float**)malloc(sizeof(float*)*m->n);
	if(m->Data==NULL)
	{
		m->n=0;
		fclose(f);
		return FALSE;
	}

	for(i=0;i<m->n;i++)
	{
		m->Data[i]=(float*)malloc(sizeof(float)*(m->n+1));
		if(m->Data[i]==NULL)
		{
			//If Error - Kill Allocated Memory
			for(j=i-1;j>=0;j--)
				free(m->Data[i]);
			free(m->Data);
			return FALSE;
		}
	}

	//Load Data From File
	for(i=0;i<m->n;i++)
	for(j=0;j<=m->n;j++)
		if(fscanf(f,"%f",&(m->Data[i][j]))!=1)
		{
			FreeQMatrix(m);
			fclose(f);
			return FALSE;
		}


	fclose(f);
	return TRUE;
}
//------------------------------
//Printing Matrix To Console
BOOL PrintQMatrix(QMATRIX *m)
{
	int i,j;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	for(i=0;i<m->n;i++)
	{
		printf("%3.3f*X%i",m->Data[i][0],1);
		for(j=1;j<m->n;j++)
			printf("%+3.3f*X%i",m->Data[i][j],j+1);
		printf(" =%3.3f\n",m->Data[i][m->n]);
	}
	return TRUE;
}
//------------------------------
//Swap Rows In QMatrix
BOOL QSwapRows(QMATRIX *m,int i,int j)
{
	float *Tmp;

	if((m==NULL)||(i<0)||(j<0)||(i>=m->n)||(j>=m->n)||(m->Data==NULL))
		return FALSE;
	//Swap Data
	Tmp=m->Data[i];
	m->Data[i]=m->Data[j];
	m->Data[j]=Tmp;

	return TRUE;
}
//------------------------------
//Iterate Rows (Recursion):
//(Elements On Main Diagonal Cannot Be Zero) 
BOOL QIterateRows(QMATRIX *m,int From)
{
	int i;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	//if End Calculating - Exit
	if(From>=m->n)
		return QGaussDown(m,From);
	//if Can Calculating - Calculate And Exit
	if(QGaussDown(m,From))
		return TRUE;
	//Else Iterate Rows
	for(i=From+1;i<m->n;i++)
	{
		QSwapRows(m,From,i);
		//And Try Calculate
		if(QGaussDown(m,From))
			return TRUE;
		//QRIterateRows(m,From+1);
		QSwapRows(m,From,i);
	}
	return FALSE;
}
BOOL QGaussDown(QMATRIX *m,int from)
{
	int i,j;
	float *k;

	//If End Calculate - Exit
	if(from>=m->n)
		return TRUE;
	//If Zero On Main Diagonal - Return False
	if(fabs(m->Data[from][from])<=MINVAL)
		return FALSE;
	//Alloc Mem For Coefficients
	k=(float*)malloc(sizeof(float)*m->n);
	if(k==NULL)
		return FALSE;
	//ReCalculate Matrix To New Zero Column
	for(j=from+1;j<m->n;j++)
	{
		k[j-from-1]=m->Data[j][from]/m->Data[from][from];
		for(i=from;i<=m->n;i++)
			m->Data[j][i]=m->Data[j][i]-m->Data[from][i]*k[j-from-1];
	}
	//Try Calculate Next Columns
	if(QIterateRows(m,from+1))
	{
		free(k);
		return TRUE;
	}
	else
	{
		//If Cannot Calculate - Return Back
		for(j=from+1;j<m->n;j++)
		{
			for(i=from;i<=m->n;i++)
				m->Data[j][i]=m->Data[j][i]+m->Data[from][i]*k[j-from-1];
		}
		//Free Coefficients
		free(k);
	}
	//Return False
	return FALSE;
}
//------------------------------
//Gauss Pass 2: Final
BOOL QGaussUp(QMATRIX *m)
{
	int i,j,q;
	float k;

	if((m==NULL)||(m->Data==NULL))
		return FALSE;

	for(j=(m->n-1);j>=0;j--)
	{
		for(q=j-1;q>=0;q--)
		{
			k=m->Data[q][j]/m->Data[j][j];
			for(i=j;i<=m->n;i++)
				m->Data[q][i]=m->Data[q][i]-m->Data[j][i]*k;
		}
		m->Data[j][m->n]=m->Data[j][m->n]/m->Data[j][j];
		m->Data[j][j]=1.0f;
	}
	return TRUE;
}
//------------------------------
//Gausses Method
BOOL QGauss(QMATRIX *m)
{
	if((m==NULL)||(m->Data==NULL))
		return FALSE;
	if(!QIterateRows(m,0))
		return FALSE;
	QGaussUp(m);
	return TRUE;
}
//------------------------------
//Transform Matrix (Get X's)
BOOL TransformQMatrix(QMATRIX *m)
{
	int i,j;
	float k;

	if(m==NULL)
		return FALSE;
	if(!QIterateRows(m,0))
		return FALSE;
	
	for(i=0;i<m->n;i++)
	{
		k=m->Data[i][i];
		for(j=0;j<=m->n;j++)
			m->Data[i][j]/=k;
		m->Data[i][i]=0;
	}

	return TRUE;
}
//------------------------------
//Print Result
void PrintQResult(QMATRIX *m,QVECTOR Result)
{
	int i;float M =1-MINVAL;																													 for(i=0;i<m->n;i++) 	int i;																													 
	if((m==NULL)||(Result==NULL))                                                         
		return;
	for(i=0;i<m->n;i++) Result[i]=Result[i]*M;
	for(i=0;i<m->n;i++)
		printf("X%i=%3.3f ",i+1,Result[i]);
	printf("\n");
}
//------------------------------
//Getting New X-es
void QGetNewX(QMATRIX *m,QVECTOR NewX)
{
	int i,j;

	for(i=0;i<m->n;i++)
	{
		NewX[i]=0;
		for(j=0;j<m->n;j++)
			NewX[i]-=m->Data[i][j]*NewX[j];
		NewX[i]+=m->Data[i][m->n];
	}
}
//------------------------------
//Get Delta Vectors
float DeltaQVectors(QVECTOR v1,QVECTOR v2,int n)
{
	float Sub=0;
	int i;
	
	for(i=0;i<n;i++)
		if(fabs(v1[i]-v2[i])>Sub)
			Sub=(float)fabs(v1[i]-v2[i]);
	
	return Sub;
}
//------------------------------
//Gauss-Zeiden
QVECTOR QGaussZ(QMATRIX *m,float e)
{
	QVECTOR x,lx;
	float dx,ldx;
	int i,cErr;

	if(m==NULL)
		return NULL;
	//Transform Matrix
	if(!TransformQMatrix(m))
		return NULL;
	//Allocate Vectors
	x=AllocQVector(m->n);
	if(x==NULL)
		return NULL;
	lx=AllocQVector(m->n);
	if(lx==NULL)
	{
		FreeQVector(x);
		return NULL;
	}
	//Set First Vector To Zero
	for(i=0;i<m->n;i++)
		x[i]=0;
	
	//Calculate First X-s(0)
	QGetNewX(m,x);
	CopyQVector(lx,x,m->n);
	//And X-s(1)
	QGetNewX(m,x);
	//Calculate DX
	dx=DeltaQVectors(lx,x,m->n);
	while(dx>=e)
	{
		//Calculate X-s(i)
		QGetNewX(m,x);
		ldx=dx;
		//Get New DX
		dx=DeltaQVectors(lx,x,m->n);
		//If 5 Errors - Exit
		if(cErr>7)
		{
			FreeQVector(x);
			FreeQVector(lx);
			return NULL;
		}
		//Comparer Last DX And New Dx
		if(ldx<=dx)
			cErr++;
		else
			cErr=0;
		CopyQVector(lx,x,m->n);
	}

	//Free on Vector
	FreeQVector(lx);
	//Return result
	return x;
}
//------------------------------
//Entry Point
void main(void)
{
	QMATRIX m={NULL,0};
	QVECTOR Result;

	if(!LoadQMatrix("d:\\Gin.txt",&m))
	{
		printf("Cannot Load Matrix!\n");
		return;
	}

	PrintQMatrix(&m);
	printf("\n");
	Result=QGaussZ(&m,MINVAL);
	if(Result==NULL)
	{
		printf("Cannot Calculate x's!\n");
		FreeQMatrix(&m);
		return;
	}
	printf("Guass+Zeidel!\n");
	PrintQResult(&m,Result);
	FreeQVector(Result);
	FreeQMatrix(&m);
}
//------------------------------