// CMPSC 201
// Gaussian Elimination Demo #1
// http://www-users.itlabs.umn.edu/classes/Summer-2002/csci1113/gauss_elim.html

#include <iostream>
#include <iomanip>
using namespace std;

const int COLS = 20;
const int ROWS = 20;

void GaussElim(double[][COLS], double[], int, int);

int main()
{
	double augment[ROWS][COLS] = {0};
	double x[ROWS] = {0};
	int i = 0;
	int j = 0;
	int rows = 0;
	int cols = 0;

	cout << "Enter number of rows and columns in augmented matrix: ";
	cin >> rows;
	cin >> cols;
	cout << "Enter augmented matrix rowwise:" << endl;
	
	for(i = 0; i < rows; i++)
	{

		for(j = 0; j < cols; j++)
		{
			cin >> augment[i][j];
		}

	}

	cout << endl << "Augmented matrix:" << endl;
	
	for(i = 0; i < rows; i++)
	{
		for(j = 0; j < cols; j++)
		{
			cout << setw(6) << augment[i][j];
		}

		cout << endl;
	}

	GaussElim(augment, x, rows, cols);
	
	cout << endl << "The linear system solution is:" << endl;
	
	for(i = 0; i < rows; i++)
	{
		cout << "x[" << i << "] = " << x[i] << endl;
	}

	return 0;
}// end of main

void GaussElim( double M[][COLS], double x[] , int rows, int cols)
{
	int i = 0;
	int j = 0;
	int k = 0;
	double multiplier = 0.0;
	
	// forward elimination
	for (i = 0; i < rows-1; i++)
	{
		for (j = i + 1; j < rows; j++)
		{
			multiplier = M[j][i]/M[i][i];

			for(k = i; k < cols ; k++)
			{
				M[j][k] = M[j][k] - multiplier * M[i][k];
			}

		}

	}

	cout << endl << "The matrix before back substitution is:" << endl;
	
	for (i = 0; i < rows; i++)
	{
		
		for (j = 0; j < cols; j++)
		{
			cout << setw(6) << M[i][j];
		}
	
		cout << endl;
	}

	// back substitution
	double sum = 0;
	
	for (int r = rows - 1; r >= 0; r--)
	{
		sum = 0;
	
		for (int c = r + 1; c < cols - 1; c++)
		{
			sum += M[r][c] * x[c];
		}
	
		x[r] = 1/M[r][r] * (M[r][cols - 1] - sum);
	}
}// end of GaussElim