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

using namespace std;

void showmatrix();
void zeilen_vertauschen(int zeile1, int zeile2);
float fround(float zahl, float stellen);

float** lgs;
float* tempzeile;
float* loesungsvektor;
int zeilen=0;

void zeilen_vertauschen(int zeile1, int zeile2)
{
	tempzeile = lgs[zeile1];
	lgs[zeile1] = lgs[zeile2];
	lgs[zeile2] = tempzeile;
}

float fround(float zahl, float stellen)
{
	return(round(zahl * pow(10,stellen) )*pow(10,-stellen));
}

void showmatrix()
{
	for(int i=0;i< zeilen;i++)
	{
		for(int j=0;j< zeilen;j++)
		{
			cout  << setw(9) << setprecision(3)<< lgs[i][j] << " ";
		}
		cout << "| " << setprecision(6) << lgs[i][zeilen] << "\n";
	}
}

int main()
{
	
	string dateiname = "";
	bool spaltenpivot = true;
	bool einzel = false;
	ifstream quelldatei;
	char puffer[255] = "";
	
	cout << "Aufgabe 14: Gleichungssysteme - Gauss-Algorithmus\n";
	cout << "Ok, welche Datei soll denn eingelesen werden? [a14-lgs1.dat]: ";
	cin.getline(puffer,sizeof(puffer),'\n');
	dateiname = (string)puffer;
	if(dateiname == "") dateiname = "a14-lgs1.dat";
	
	// Oeffnen der Datei
	quelldatei.open(dateiname.c_str(), ios::in);
	if(!quelldatei)
	{
		cerr << "Fehler: " << dateiname << " kann nicht geoeffnet werden.\n";
		exit(1);
	}
	
	// Einlesen der Zeilenanzahl
	quelldatei >> zeilen;
	
	if(zeilen > 0)
	{
		cout << "Unsere Matrix hat " <<  zeilen << " Zeilen\n";
	}
	else
	{
		cerr << "Fehler: Die Matrix konnte nicht eingelesen werden!\n";
		exit(2);
	}
	
	cout << "Soll Spaltenpivotisierung verwendet werden? (j/n) [j]: ";
	cin.getline(puffer,sizeof(puffer),'\n');
	if((string)puffer == "n")
	{
		spaltenpivot = false;
		cout << "Gut, dann eben ohne Spaltenpivotisierung...\n";
	}
	
	cout << "Sollen Einzelschritte angezeigt werden? (j/n) [n]: ";
	cin.getline(puffer,sizeof(puffer),'\n');
	if((string)puffer == "j")
	{
		einzel = true;
		cout << "Ok, Einzelschritte werden angezeigt...\n";
	}
	
	// Initialisieren der Matrix & der Tempzeile
	tempzeile = new float[zeilen + 1];
	lgs = new float*[zeilen];
	for(int i =0; i < zeilen;i++)
	{
			lgs[i] = new float[zeilen];
			for(int j=0;j<= zeilen;j++)
			{
				lgs[i][j] = 0;
			}
	}
	
	// Initialisierung des Loesungsvektors
	loesungsvektor = new float[zeilen];
	for(int i =0; i < zeilen;i++) loesungsvektor[i] = 0;
	
	// Einlesen der Matrix aus der Datei
	// Elemente von A
	for(int i=0; i < zeilen;i++)
	{
		for(int j=0; j < zeilen;j++)
		{
			quelldatei >> lgs[i][j];
		}
	}
	// Elemente von b
	for(int i=0;i< zeilen;i++)
	{
		quelldatei >> lgs[i][zeilen];
	}
	
	
	// Ausgabe der Matrix
	cout << "Orginalmatrix:\n";
	showmatrix();
	
	//Anwenden des Gauss-Verfahrens
	float fakt=1; // Linearer Fakator zum skalieren der Zeile
	for(int i=0;i< zeilen-1;i++)
	{
		// Bearbeiten der Spalte i:
		
		// Fuer Spaltenpivotisierung: vertauche Zeilen ,so dass element mit groesstem
		// |[j][i]| an oberster Stelle steht.
		if(spaltenpivot)
		{
			// Finde betragsgroesstes Element:
			int groesstes=0;
			float groessterwert = 0;
			for(int l = i;l < zeilen;l++)
			{
				if(groessterwert < fabs(lgs[l][i]))
				{
					groessterwert = fabs(lgs[l][i]);
					groesstes = l;
				}
			}
			//Setze Zeile mit groesstem element nach oben
			if(groesstes != i) zeilen_vertauschen(i,groesstes);
		}
		
		for(int j=i+1;j< zeilen ;j++)
		{
			// Nimm Zeile [i] und addiere sie so zu Zeile j,
			// so dass [j][i]=0
			fakt = - lgs[j][i] / lgs[i][i];
			for(int k=0;k <= zeilen;k++) 
			{
				if(k <= i)
				{
					lgs[j][k] = 0;
				}else{
					lgs[j][k] = lgs[j][k] + (fakt * lgs[i][k]);
				}
			}
		}
		//Ausgabe der Einzelschritte wenn erwuenscht
		if(einzel)
		{
			cout << "Schritt:"<<i+1<<"\n";
			showmatrix();
		}
		
	}
	
	// Ausgabe der dreiecks Matrix
	cout << "Dreieckmatrix:\n";
	showmatrix();
	
	
	// Rekursives Lösen des LGS
	float tempzahl;
	for(int i=(zeilen -1);i>= 0;i--)
	{
		tempzahl = lgs[i][zeilen];
		for(int j = (zeilen -1); j > i; j--)
		{
			tempzahl -= lgs[i][j] * loesungsvektor[j];
		}
		loesungsvektor[i] = tempzahl / lgs[i][i];
	}
	
	// Ausgabe des Lösungsvektors
	cout << "Der Lösungsvektor lautet (auf 2 Nachkommastellen gerundet):\n";
	for(int i=0; i< zeilen;i++)
	{
		cout << "x" << i+1 << "= " << fround(loesungsvektor[i],2) << "\n";
	}
}