/************************************************/
/* ECE-1011 Spring 2003                         */
/* Section:......... 0                          */
/* Assignment:...... 8                          */
/* Programmer:...... BAHN, William L.           */
/* File Contents:... SOURCE CODE                */
/* Filename:........ s80_key                    */
/* Due Date:........ 31 MAR 03                  */
/* E-mail Address:.. wbahn@eas.uccs.edu         */
/************************************************/

/*****************************************************************************/
/*************************** GAUSSIAN ELIMINATOR *****************************/
/*****************************************************************************/

/* TAKEN DIRECTLY FROM HOMEWORK ASSIGNMENT:

TASKS TO BE COMPLETED

1)  Ask the user for a data file name (they supply just the name and you
	 add a *.dat extension to it).
2)  Open the data file with all of the associated error checking.
3)  Determine the number of equations in the file.
4)  Report the number of equations found to the user and to the output file.
5)  Ask the user if they wish to use row-pivoting (see Prob 46 on page 299).
6)  Ask the user if they wish to use column-pivoting (see Prob 48 on page 299).
7)  Read the data from the file into an appropriate array.
8)  Solve the system of equations.
9)  Output the results.

DATA FILE FORMAT

	M Equations:
	* M lines of data
	* (M+1) floating point values on each line.
	* No more than 20 equations.
	* There is no header or trailer record.

OUTPUT FILE FORMAT

	NAME: LAST NAME, First Name  (Last Name ALL CAPS, first name normal caps)
	E-MAIL: e-mail address
	VERSION: Version#.Revision#

	DATA SET: setname
	EQUATIONS: <number of equations in data set>
	TOTAL DATA FOUND: <number of values read from file>

	OPTIONS:
	ROW PIVOTING: <Y or N> (should reflect what was used, not what was requested)
	COL PIVOTING: <Y or N> (should reflect what was used, not what was requested)

	SOLUTION:
	x0 = <value>
	x1 = <value>
	...
	x(M-1) = <value>

*/

/*****************************************************************************/
/* BASIC APPROACH                                                            */
/*****************************************************************************/
/*

The goal is to take a matrix of the form:

	a11 a12 a13 a14
	a21 a22 a23 a24
	a31 a32 a33 a34

where, for instance, the first line represents the equation:

	a11*x1 + a12*x2 + a13*x3 = a14

and solve the system of N simultaneous equations.

The goal is therefore to manipulate the matrix into the final form of:

	 1 0 0 sx1
	 0 1 0 sx2
	 0 0 1 sx3

Based on working a hand example, the following tasks must be doable:

1) Normalize an equation to a specified column number - in other words,
	force the coefficient in Column C to be equal to 1.0. Note that if the
	equation contains a zero in that column it can't be nomalized.

2) Given two equations normalized to the same column, eliminate that
	column from one of the equations. Note that if one of the equations
	contains a zero in that column it is a special case.

3) Back substitute the value for one variable to eliminate it from the
	remaining equations.

Psuedocode to Normalize an equation:

	To normalize equation R to column C:

		1) If the coefficient in column C is zero, can't be normalized.
		2) For each column in the row R:
			2.1) Divide the coefficient by the coefficient in column C.

Pseudocode to Eliminate a variable from two equations:

	To eliminate variable in column C from equation R using:

		1) If the coefficient in row R column C is zero, it is already done.
		2) Normalize Row C to Column C.
		3) Normalize Row R to Column C.

		We now have:
					  C
				a11 a12 a13 a14
				a21  1  a23 a24
			R	a31  1  a33 a34

		4) For each column c in row R :
			4.1) data[R[[c] -= data[C][c]

Pseudocode to Back Substitute a variable:

	To Back Substitute column C to eliminate it from row R:

		We have:
							C
			R  a11 a12 a13 a14
				a21  1  a23 a24
				 0   0   1  s34

		Since row R represents the equation:

			a11*x1 + a12*x2 + a13*x3 = a14

		And since row C informs us that:

			x3 = s34

		We have:

			a11*x1 + a12*x2 +  0  =  a14 - a13*s34

		Hence:

		1) data[R][N] -= data[R][C] * data[C][N];
		2) data[R][C] = 0;

*/

/*============== STANDARD LIBRARY INCLUDE FILES =============================*/

#include <stdio.h>  /* printf(), scanf(), fprintf(), fscanf() */
#include <string.h> /* strcpy(), strcat() */
#include <math.h>   /* fabs() */

/*============== PROGRAMMER DEFINED CONSTANTS ===============================*/

/*****************************************************************************/
/* Version and Author Control                                                */
/*****************************************************************************/

#define PROGRAMTITLE     ("Gaussian Eliminator")
#define VERSION          (8)
#define REVISION         (1)
#define PROGRAMMER       ("BAHN, William")
#define PROGRAMMER_EMAIL ("wbahn@eas.uccs.edu")
#define LASTMODDATE      ("29 MAR 2003")

/*****************************************************************************/
/* Programmer Defined Constants                                              */
/*****************************************************************************/

#define TRUE  (1==1)
#define FALSE (!TRUE)
#define BLANKLINE printf("\n")
#define MAXEQUATIONS (20)
#define COLS (MAXEQUATIONS+1)
#define EPSILON (0.000001)

/* EXIT CODES */
#define BAD_DATAFILE   (1)
#define BAD_OUTPUTFILE (2)
#define BAD_EQUATIONS  (3)

/*****************************************************************************/
/* Function Prototypes                                                       */
/*****************************************************************************/

void  PrintHeader(void);
void  PrintOutputFileHeader(FILE *output_fp, char dataset[]);
void  GetDataSetName(char dataset[]);
FILE *OpenDataFile(char dataset[]);
FILE *OpenOutputFile(char dataset[]);
int   CountEquations(FILE *infp, FILE *outfp);
int   GetEquationsFromDataFile(FILE *fp, double d[][COLS], int o[], int eqns);
void  PrintOptions(FILE *fp, int rflag, int cflag);
int   PrintEquations(double d[][COLS], int order[], int eqns);
int   DoRowPivot(void);
int   DoColPivot(void);
int   SolveEqns(double d[][COLS], int order[], int eqns, int rpiv, int cpiv);
int   PivotRow(double d[][COLS], int row, int eqns);
int   PivotCol(double d[][COLS], int order[], int col, int eqns);
int   NormalizeRowToCol(double d[][COLS], int row, int col, int eqns);
int   EliminateColFromRow(double d[][COLS], int row, int col, int eqns);
int   BackSubstituteColIntoRow(double d[][COLS], int row, int col, int eqns);
int   OutputSolutions(FILE *fp, double d[][COLS], int order[], double eqns);
int   getresponse(char msg[], char s[]);

/*****************************************************************************/
/* MAIN FUNCTION                                                             */
/*****************************************************************************/

int main(void)
{
	char dataset[128];
	double data[MAXEQUATIONS][COLS];
	int order[MAXEQUATIONS];
	int eqns;
	int rowpivflag, colpivflag;
	FILE *data_fp, *output_fp;

	PrintHeader();

	GetDataSetName(dataset);                                 /* TASK  #1    */

	if(NULL == (data_fp = OpenDataFile(dataset)))            /* TASK  #2    */
	{
		fcloseall();
		return(BAD_DATAFILE);
	}

	if(NULL == (output_fp = OpenOutputFile(dataset)))        /* TASK  #2    */
	{
		fcloseall();
		return(BAD_OUTPUTFILE);
	}

	PrintOutputFileHeader(output_fp, dataset);

	eqns = CountEquations(data_fp, output_fp);               /* TASKS #3,#4 */

	if(!eqns)
	{
		fcloseall();
		return(BAD_EQUATIONS);
	}

	printf("Resetting Data file - Closing and reopening.\n");
	fclose(data_fp);                                         /* Reset File  */
	if(NULL == (data_fp = OpenDataFile(dataset)))
	{
		fcloseall();
		return(BAD_DATAFILE);
	}

	rowpivflag = DoRowPivot();                               /* TASK  #5    */
	colpivflag = DoColPivot();                               /* TASK  #6    */
	PrintOptions(output_fp, rowpivflag, colpivflag);
	GetEquationsFromDataFile(data_fp, data, order, eqns);    /* TASK  #7    */
	PrintEquations(data, order, eqns);
	if(0==SolveEqns(data, order, eqns, rowpivflag, colpivflag)) /* TASK  #8  */
	{
		PrintEquations(data, order, eqns);
		OutputSolutions(output_fp, data, order, eqns);        /* TASK  #9    */
	}
	else
	{
		printf (           "No Solution Found!\n");
		fprintf(output_fp, "No Solution Found!\n");
	}

	fcloseall();
	return(0);
}

/*****************************************************************************/
/* SUPPORT FUNCTIONS                                                         */
/*****************************************************************************/

void PrintHeader(void)
{
	printf("===============================================================\n");
	printf("%s Version %i.%i\n", PROGRAMTITLE, VERSION, REVISION);
	printf("Programmer: %s\n", PROGRAMMER);
	printf("Date of Last Mod: %s\n", LASTMODDATE);
	printf("===============================================================\n");
	BLANKLINE;
	return;
}

void GetDataSetName(char dataset[])
{

	printf("Enter the name of the data set (8 characters max): ");
	scanf("%s", dataset);

	dataset[8] = '\0';   /* Force the data set name to end at 8 characters */
	printf("Data Set: %s\n", dataset);
	BLANKLINE;

	return;
}

FILE *OpenDataFile(char dataset[])
{
	char filename[16];
	FILE *fp;

	/* Create the name for the input data file, open and verify */
	strcpy(filename, dataset);
	strcat(filename, ".dat");

	if( NULL == (fp = fopen(filename, "rt")) )
	{
		printf("Failure to open data file <%s> for reading.\n", filename);
		printf("No attempt has or will be made to open output file.\n");
		printf("ABORTING PROGRAM!\n");
		BLANKLINE;
	}
	else
	{
		printf("Data file <%s> successfully opened for reading.\n", filename);
		BLANKLINE;
	}

	return(fp);
}

FILE *OpenOutputFile(char dataset[])
{
	char filename[16];
	FILE *fp;

	/* Create the name for the output file, open and verify */

	strcpy(filename, dataset);
	strcat(filename, ".out");

	if( NULL == 	(fp = fopen(filename, "wt")) )
	{
		printf("Failure to open output file <%s> for writing\n", filename);
		printf("ABORTING PROGRAM!\n");
		BLANKLINE;
	}
	else
	{
		printf("Output File <%s> successfully opened for writing.\n", filename);
		BLANKLINE;
	}

	return(fp);
}

void PrintOutputFileHeader(FILE *output_fp, char dataset[])
{
	BLANKLINE;
	fprintf(output_fp, "NAME: %s\n", PROGRAMMER);
	fprintf(output_fp, "E-MAIL: %s\n", PROGRAMMER_EMAIL);
	fprintf(output_fp, "VERSION: %i.%i\n", VERSION, REVISION);
	fprintf(output_fp, "\n");
	fprintf(output_fp, "DATA SET: %s\n", dataset);
	BLANKLINE;
	return;
}

int   CountEquations(FILE *infp, FILE *outfp)
{
	int values, eqns;
	double value;

	values = 0;
	while(1 == fscanf(infp, "%lf", &value))
		values++;

	eqns = MAXEQUATIONS; 
	while( ((eqns)*(eqns+1) != values) && (eqns > 0) ) 
		eqns--;

	printf (       "EQUATIONS: %i\n", eqns);
	fprintf(outfp, "EQUATIONS: %i\n", eqns);

	printf (       "TOTAL DATA FOUND: %i\n\n", values);
	fprintf(outfp, "TOTAL DATA FOUND: %i\n\n", values);

	return(eqns);
}

int GetEquationsFromDataFile(FILE *fp, double d[][COLS], int order[], int eqns)
{
	int row, col;

	for(row = 0; row < eqns; row++)
	{
		order[row] = row;
		for(col = 0; col <= eqns; col++)
			fscanf(fp, "%lf", &d[row][col]);
	}

	return(0);
}

int PrintEquations(double d[][COLS], int order[], int eqns)
{
	int row, col;

	for(col = 0; col < eqns; col++)
		printf("   %2i  ", order[col]);
	printf("   =   \n");

	for(row = 0; row < eqns; row++)
	{
		for(col = 0; col <= eqns; col++)
			printf("%7.2f", d[row][col]);
		printf("\n");
	}
	BLANKLINE;

	return(0);
}

int   DoRowPivot(void)
{
	int res;

	res = getresponse("Do you want to perform row pivoting? ", "yYnN");

	switch(res)
	{
		case 'y':
		case 'Y': return(TRUE);

		case 'n':
		case 'N': return(FALSE);
	}

	return(FALSE); /* Should never reach this point */
}

int   DoColPivot(void)
{
	int res;

	res = getresponse("Do you want to perform column pivoting? ", "yYnN");

	switch(res)
	{
		case 'y':
		case 'Y': return(TRUE);

		case 'n':
		case 'N': return(FALSE);
	}

	return(FALSE); /* Should never reach this point */
}

void PrintOptions(FILE *fp, int rflag, int cflag)
{
	printf (    "OPTIONS:\n");
	fprintf(fp, "OPTIONS:\n");

	printf (    "ROW PIVOTING: ");
	fprintf(fp, "ROW PIVOTING: ");

	if(rflag)
	{
		printf(     "Y\n");
		fprintf(fp, "Y\n");
	}
	else
	{
		printf(     "N\n");
		fprintf(fp, "N\n");
	}

	printf (    "COL PIVOTING: ");
	fprintf(fp, "COL PIVOTING: ");

	if(cflag)
	{
		printf(     "Y\n");
		fprintf(fp, "Y\n");
	}
	else
	{
		printf(     "N\n");
		fprintf(fp, "N\n");
	}

	printf(     "\n");
	fprintf(fp, "\n");

	return;
}


int   SolveEqns(double d[][COLS], int order[], int eqns, int rpiv, int cpiv)
{
	int r, c;
	int rowpivots, colpivots;

	rowpivots = 0;
	colpivots = 0;

	/* Eliminate Lower Left Half of Matrix */
	for(c = 0; c < eqns; c++)
	{
		if(rpiv)
			rowpivots += PivotRow(d, c, eqns);

		if(cpiv)
			colpivots += PivotCol(d, order, c, eqns);

		if(0 != NormalizeRowToCol(d, c, c, eqns))
			return(-1);
		for(r = c+1; r < eqns; r++)
		{
			if(0 != NormalizeRowToCol(d, r, c, eqns))
				return(-2);
			EliminateColFromRow(d, r, c, eqns);
		}
	}

	/* Back Substitute Upper Right Half of Matrix */
	for(c = eqns - 1; c > 0; c--)
		for(r = c-1; r >= 0; r--)
			BackSubstituteColIntoRow(d, r, c, eqns);

	return(0);
}

int PivotRow(double d[][COLS], int row, int eqns)
{
	/*	To perform Row Pivot about Row R:

		1) Find the largest (absolute value) coefficient in Col R from Row R on.
		2) Swap row R with the row found in Step 1.

	*/

	int c, r, maxrow;
	double temp;

	maxrow = row;
	for(r = row+1; r < eqns; r++)
	{
		if(fabs(d[r][row]) > fabs(d[maxrow][row]))
			maxrow = r;
	}

	if(maxrow != row)
	{
		for(c = 0; c <= eqns; c++)
		{
			temp = d[row][c];
			d[row][c] = d[maxrow][c];
			d[maxrow][c] = temp;
		}
		return(1);
	}
	return(0);
}

int PivotCol(double d[][COLS], int order[], int col, int eqns)
{
	/*	To perform Col Pivot about Col C:

		1) Find the largest (absolute value) coefficient in Row C from Col C on.
		2) Swap col R with the col found in Step 1.
		3) Keep track of the order of the variables.

			The value stored in order[i] is the column in d[] where variable i is.
	*/

	int c, r, maxcol;
	double temp;

	maxcol = col;
	for(c = col+1; c < eqns; c++)
	{
		if(fabs(d[col][c]) > fabs(d[col][maxcol]))
			maxcol = c;
	}

	if(maxcol != col)
	{
		for(r = 0; r < eqns; r++)
		{
			/* Swap columns */
			temp = d[r][col];
			d[r][col] = d[r][maxcol];
			d[r][maxcol] = temp;
			/* Swap variable maps */
		}
		order[order[col]] = maxcol;
		order[order[maxcol]] = col;
		return(1);
	}
	return(0);
}

int NormalizeRowToCol(double d[][COLS], int row, int col, int eqns)
{
	/* To normalize equation R to column C:

		1) If the coefficient in column C is zero, can't be normalized.
		2) For each column in the row R:
			2.1) Divide the coefficient by the coefficient in column C.

	*/

	int c;
	double norm;

	if(fabs(d[row][col]) <= EPSILON)
		return(-1);

	for(c = 0, norm = d[row][col]; c <= eqns; c++)
		d[row][c] /= norm;

	return(0);
}

int EliminateColFromRow(double d[][COLS], int row, int col, int eqns)
{
	/*	To eliminate variable in column C from equation R:

		1) If the coefficient in row R column C is zero, it is already done.
		2) Normalize Row C to Column C. - ASSUME DONE AT HIGHER LEVEL
		3) Normalize Row R to Column C. - ASSUME DONE AT HIGHER LEVEL
		4) For each column c in row R :
			4.1) data[R[[c] -= data[C][c]
	*/

	int c;

	if(fabs(d[row][col]) <= EPSILON)
		return(0);

	for(c = 0; c <= eqns; c++)
		d[row][c] -= d[col][c];

	return(0);
}

int BackSubstituteColIntoRow(double d[][COLS], int row, int col, int eqns)
{

	/* To Back Substitute column C to eliminate it from row R:

		We have:
							C
			R  a11 a12 a13 a14
				a21  1  a23 a24
				 0   0   1  s34

		Since row R represents the equation:

			a11*x1 + a12*x2 + a13*x3 = a14

		And since row C informs us that:

			x3 = s34

		We have:

			a11*x1 + a12*x2 +  0  =  a14 - a13*s34

		Hence:

		1) data[R][N] -= data[R][C] * data[C][N];
		2) data[R][C] = 0;

	*/

	d[row][eqns] -= d[row][col] * d[col][eqns];
	d[row][col] = 0.0;

	return(0);
}

int   OutputSolutions(FILE *fp, double d[][COLS], int order[], double eqns)
{
	int eq;

	for(eq = 0; eq < eqns; eq++)
	{
		printf (    "x%i = %f\n", eq, d[order[eq]][eqns]);
		fprintf(fp, "x%i = %f\n", eq, d[order[eq]][eqns]);
	}

	return(0);
}

int getresponse(char msg[], char s[])
{
	int res, i, max;
	int exitflag;

	exitflag = FALSE;
	max = strlen(s);
	while(!exitflag)
	{
		printf("%s", msg);
		while('\n' == (res = getchar())); /* Get 1st non-CR from buffer */
		i = 0;
		while( (i < max) && (res != s[i++]) ); /* Look for 'res' in s[] */
		if(i < max)
			exitflag = TRUE;                    /* 'res' found in s[]    */
		else
			printf("Invalid Response! - Please Reenter\n");
		while('\n' != getchar());                /* Clear Buffer regardless */
	}

	return(res);
}