#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cmath>

using namespace std;

void rysujMacierz(double ** elem, int dim, int dim2, int first_col = 0)
{
	cout << setprecision(7) << fixed;

	for(int i = 0; i < dim; i++)
	{
		for(int j = first_col; j < dim2; j++)
		{
			cout << setw(12) << elem[i][j] << " ";
		}
		cout << endl;
	}
}

void rysujWektor(double * elem, int dim)
{
	cout << setprecision(7) << fixed;

	for(int i = 0; i < dim; i++)
	{
		cout << setw(12) << elem[i] << " ";
	}
	cout << endl;
}

int main()
{
	int dim;
	double ** macierz; //glowna macierz
	double * wektorB; //wektor wyrazow wolnych
	double * x; //rozwiazanie
	double * y;

	dim = 4;

	macierz = new double * [dim];
	wektorB = new double [dim];
	x = new double [dim + 1];
	y = new double [dim + 1];

	for(int i = 0; i < dim; i++)
	{
		macierz[i] = new double [2 * dim];
	}

	macierz[0][4] = 1.0;
	macierz[0][5] = 20.0;
	macierz[0][6] = -30.0;
	macierz[0][7] = -4.0;

	macierz[1][4] = 4.0;
	macierz[1][5] = 20.0;
	macierz[1][6] = -6.0;
	macierz[1][7] = 50.0;

	macierz[2][4] = 9.0;
	macierz[2][5] = -18.0;
	macierz[2][6] = 12.0;
	macierz[2][7] = -11.0;

	macierz[3][4] = 16.0;
	macierz[3][5] = -15.0;
	macierz[3][6] = 14.0;
	macierz[3][7] = 130.0;

	wektorB[0] = 0.0;
	wektorB[1] = 114.0;
	wektorB[2] = -5.0;
	wektorB[3] = 177.0;

	cout << "\npodana macierz wspolczynnikow:" << endl;

	rysujMacierz(macierz, dim, 2 * dim, dim);

	cout << "\npodany wektor wyrazow wolnych:" << endl;

	rysujWektor(wektorB, dim);


	//dolaczamy macierz jednostkowa
	for(int i = 0; i < dim; i++)
	{
		for(int j = 0; j < dim; j++)
		{
			if (i == j)
			{
				macierz[i][j] = 1.0;
			}
			else
			{
				macierz[i][j] = 0.0;
			}
		}
	}

	cout << "\nMacierz klatkowa:" << endl;

	rysujMacierz(macierz, dim, 2 * dim);


	//rozklad LU
	double e = 1e-9; //0
	double m;

	int i;
	int j;
	int k;

	for(i = 0; i < dim; i++)
	{

		//czesciowy wybor elementu podstawowego
		if(fabsl(macierz[i][dim + i]) < e)
		{
			cout << "element na przekatnej zbyt zblizony zera!" << endl;

			double wektor_mod_max = fabsl(macierz[i][dim + i]); // element o maksymalnym module.
			int k_max = i; // numer wiersza w ktorym znaleziono maks element

			// szukamy elementu o maks. module
			// w elementach z kolumny i.
			for(int k = i; k < dim; k++)
			{
				if(fabsl(macierz[k][dim + i]) > wektor_mod_max)
				{
					wektor_mod_max = fabsl(macierz[k][dim + i]);
					k_max = k;
				}
			}

			// jesli znaleziony element ma modul 0, konczymy.
			if(wektor_mod_max < e)
			{
				cout << "\nnie znaleziono elementu podstawowego roznego od zera - rozw. nie istnieje!";
				exit(-1);
			}

			// jesli znaleziono element o module wiekszym niz 0, zamieniamy wiersze.
			// zamieniamy k_max-ty wiersz z i-tym wierszem.

			double tmp;

			for(int k = 0; k < i; k++)
			{
				// zamieniamy tylko fragmenty wiersza mac. L.
				tmp = macierz[k_max][k];
				macierz[k_max][k] = macierz[i][k];
				macierz[i][k] = tmp;
			}

			for(int k = dim; k < 2 * dim; k++)
			{
				// Zamianiamy wiersz mac. U.
				tmp = macierz[k_max][k];
				macierz[k_max][k] = macierz[i][k];
				macierz[i][k] = tmp;
			}

			cout << "zamieniono elementy z wiersza " << k_max + 1 << " z wierszem " << i + 1 << "." << endl;

			rysujMacierz(macierz, dim, 2 * dim);
		}


		//gauss LU
		for(j = i + 1; j < dim; j++)
		{
			m = macierz[j][dim + i] / macierz[i][dim + i];

			cout << "\nwspolczynnik m = " << macierz[j][dim + i] << " / " << macierz[i][dim + i] << "." << endl;

			for(k = dim; k >= 0; k--)
			{
				macierz[j][dim + k] -= m * macierz[i][dim + k];

				cout << "Element " << j << " " << dim + k << " -= " << m << " * " << macierz[i][dim + k] << "." << endl;
			}

			// do macierzy jednostkowej z lewej wpisujemy wspolczynniki
			// uzywane przy eliminacji Gaussa.

			macierz[j][i] = m;
		}

		cout << "\nMacierz po kroku:" << endl;
		rysujMacierz(macierz, dim, 2 * dim);
	}


	//wyswietlamy macierze L i U
	cout << "\nMacierz L:" << endl;

	rysujMacierz(macierz, dim, dim);

	cout << "\nMacierz U:" << endl;

	rysujMacierz(macierz, dim, 2 * dim, dim);

	//obliczanie y i x

	//y
	y[0] = wektorB[0];

	for(i = 1; i < dim; i++)
	{
		y[i] = wektorB[i];
		for(j = 1; j < i + 1; j++)
		{
			y[i] = y[i] - macierz[i][j-1] * y[j - 1];
		}
	}

	cout << "wektor y:" << endl;

	rysujWektor(y, dim);

	//x
	x[dim - 1] = y[dim - 1] / macierz[dim - 1][dim - 1 + dim];

	for(i = dim - 1; i >= 0; i--)
	{
		x[i] = y[i];

		for(j = i + 1; j < dim; j++)
		{
			x[i] = x[i] - macierz[i][j + dim] * x[j];
		}
		x[i] = x[i] / macierz[i][i + dim];
	}

	cout << "Wektor x:" << endl;

	rysujWektor(x, dim);


	system("PAUSE");
	return 0;
	}