/*
 *
 *		Gaussian Elimination
 *
 *		Duy Nguyen
 */



#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>


#define MAX 20



/* Function prototypes */

void normalize(double matrix[MAX][MAX], double inverse[MAX][MAX], int size);
int upper_triangulate(double matrix[MAX][MAX], double inverse[MAX][MAX], 
	int size);
void Switch(double matrix[MAX][MAX], int row, int row_move, int size);
void perform_elt(double matrix[MAX][MAX], double inverse[MAX][MAX], int row, 
	int size);
void lower_triangulate(double matrix[MAX][MAX], double inverse[MAX][MAX], 
	int size);




main()
{
	double matrix[MAX][MAX], inverse[MAX][MAX];
	int i,j, singular, size;
	char filename[12];
	int count;
	ifstream ifp;
	ofstream ofp;



     /* Get input file.  File format should contain:
	integer specifying the dimension of the square matrix
	row 1 of the matrix
	row 2 of the matrix
	etc...
      */

	cout << "What is the name of the input file?  ";
	cin.getline(filename, 12);
	cout << "\n";

     /* Inform user if the file does not exist */
	ifp.open(filename, ios::in);
	if (!ifp)
	{	
	   cout << "Error opening " << filename << ".\n";
	   exit(1);
	}

     /* Get the 1st element of the file =  size of matrix */
	ifp >> size;

     /* Get the rest of the file and initialize the inverse matrix to I */
	for (i=0; i<size; i++)
	{	
	   for(j=0; j<size; j++)
	   {
		 ifp >> matrix[i][j];
	      inverse[i][j]=0;
	   }
	   inverse[i][i]=1.0;
	}
	ifp.close();

     /* Triangulate the original matrix into an upper triangle using
  	elementary row operations.  Perform the same exact operations
	on the Inverse matrix.  The return value is whether the matrix
	in singular or not.  A singular matrix is determined when the
	determinant of the matrix is zero.  In the case of the upper
	triangular matrix, the determinant is simply the product of 
	the elements on the main diagonal.  Hence, a matrix if singular
	if one of the elements on the main diagonal is zero.  This is
	all done from within the upper_triangulate function.
       */

	singular = upper_triangulate(matrix, inverse, size);

	if(singular == 1)
	{
	   cout << "Singular matrix\n.";
	   exit(1);
	}

     /* Store results into an external file */
	ofp.open("results.dat", ios::out);

     /* Normalize the elements on the main diagonal of matrix to one.
	Perform the same operations on inverse.
      */
	normalize(matrix, inverse, size);

     /* Perform elementatry operations on the upper triangular form of
	the original matrix to get the Indentity matrix starting from
	the bottom up.  Perform the same elementary operations on the
	inverse matrix.
      */

	lower_triangulate(matrix, inverse, size);
   
     /* The inverse of the original matrix is now stored in inverse.
	Print it to a file
      */
	for (i = 0;i < size;i++)
	{
   	   for (j = 0;j < size;j++)
		    ofp << inverse[i][j];
		 ofp << "\n";
	}
        ofp.close();

}  /* main */



            

/* Normalize the main diagonal of matrix. Perform the same operations on
   inverse.  */

void normalize(double matrix[MAX][MAX], double inverse[MAX][MAX], int size)
{
	int i, j;
	double factor;

	for(i=0; i<size; i++)
	{
	   /* The element on the main diagonal */
	   factor=matrix[i][i];
	   for(j=0; j<size; j++)
	   {
	      /* Divide the whole row of each matrix by the main diagonal */
	      matrix[i][j] /= factor;
	      inverse[i][j] /= factor;
	   }
	}
}  /* normalize */



/* Triangulate a square matrix into an upper triangular matrix.  Perform
   the same operations on the inverse matrix.  Return a 0 if matrix is
   singular and 1 if it is invertable.  Singularity is determined by
   whether one of the element on the main diagonal of the upper triangular
   matrix is zero.  */

int upper_triangulate(double matrix[MAX][MAX], double inverse[MAX][MAX], 
			int size)
{
	int row, column, stop, row_move, singular;


    /* Start at the top of the matrix.  row and column keep track of the
       current pivot point.  row_move moves down the matrix to find the
       first occurence of a nonzero element on the pivot column.  */

    /* Stop is true once the matrix is completely upper triangular */

	row = column = 0;
	stop = 0;
	singular = 0;
	while (!stop)
	{
	   /* row_move starts at row, if row is zero, then it goes down
	      the column to find the first nonzero and then call the
	      Switch function to exchange the rows.  */
	   row_move = row;  
	   while (matrix[row_move][column] == 0)
	   {
	      row_move++;
	      if (row_move >= size)
	      {
		 singular=1;
	         return (singular);
	      }
	   }

	   /* Switch the row and row_move.  If they are the same, then
	      the function simply does nothing. */
	   Switch(matrix, row, row_move, size);
	   Switch(inverse, row, row_move, size);

	   /* Perform elementatry operations on matrix until all the elements
	      below the pivot column are zeroed out.  Perform the same
	      operations on the inverse matrix. */
	   perform_elt(matrix, inverse, row, size);
	   row++;
	   column++;

	   /* We've completely upper triangulated matrix. */
	   if (row >= size)
	      stop = 1;
	}  /* while */

	return(singular);	

}  /* upper_triangulate */	

                                      

/* Switch two rows of a matrix */

void Switch(double matrix[MAX][MAX], int row, int row_move, int size)
{
	int i;
	double temp[3];


	if (row == row_move)
	   return;

	/* Store one row in a temporary vector first, place the other row
	   into this row, and store the temporary vector into the other row. */

	for(i = 0; i < size; i++)
	{
	   temp[i] = matrix[row_move][i];
	   matrix[row][i] = matrix[row_move][i];
	}
	for(i=0; i<size; i++)
	   matrix[row_move][i] = temp[i];

} /* Switch */ 




/* Perform elementary operations on a matrix to zero out elements below the
   pivoted column */

void perform_elt(double matrix[MAX][MAX], double inverse[MAX][MAX], int row, 
		 int size)
{
	int i, j;
	double factor;


	/* Divide the elements below the pivoted column by the pivoted 
	   element.  Multiply the pivoted row by that element times -1
	   and added to the row containing the element. */

	for (i = row+1; i < size; i++)
	{
	   factor = -1.0*matrix[i][row]/matrix[row][row];
	   for (j = row; j < size; j++)
	      matrix[i][j] += factor*matrix[row][j];
	   for (j = 0;j < size;j++)
	      inverse[i][j] += factor*inverse[row][j];
	}

} /* perform_elt */



/* Lower triangulate a matrix starting from the bottom up. */

void lower_triangulate(double matrix[MAX][MAX], double inverse[MAX][MAX], 
			int size)
{
	int i, j, k;
	double factor;


	for (i = size-1;i > 0;i--)
	   for (k = i-1;k >= 0;k--)
	   {
	      factor = -1.0*matrix[k][i]/matrix[i][i];
	      for (j = k;j >= 0;j--)
	         matrix[k][j] += factor*matrix[i][j];
	      for (j = size-1;j >= 0;j--)
		 inverse[k][j] += factor*inverse[i][j];
	   }

}