// MO LAB 5 LU.cpp : Defines the entry point for the console application.
//


#include "stdafx.h"
#include<iostream>
#include <cmath>
#include <fstream>
#include <cstdlib>
#define N 4
using namespace std;

double **stworzMacierz(int n);
void usunMacierz(double **matrix, int n);
void uzupelnijMacierz(double **matrix, int n);
void pokarzMacierz(double **matrix, int n);
int  wybor_elementu_podstawowego(double **matrix, int n, int j, int *index);
void gauss(double **matrix, int *index);
void macierzL(double**l, double *b, int *index, int n);
void macierzU(double**u, double *b, int *index, int n);
void pokarzWektor(double *b, int n);

int main()
{

	double** matrix = stworzMacierz(N);
	double b[4] = { 0.0, 114.0, -5.0, 177.0 };
	int index[4] = { 0, 1, 2, 3 };

	uzupelnijMacierz(matrix, N);
	cout << endl << "Macierz A: " << endl;
	pokarzMacierz(matrix, N);

	cout << endl  << "Wektor B: " << endl;
	pokarzWektor(b, N);

	gauss(matrix, index);
	cout << endl <<endl<< "Macierz po eliminacji Gaussa: " << endl;
	pokarzMacierz(matrix, N);

	macierzL(matrix, b, index, N-1);
	cout << endl << endl << "Wektor Y" << endl;;

	pokarzWektor(b, N);
	cout << endl << endl;
	macierzU(matrix, b, index, N-1);

	cout << "Wektor X " << endl;
	pokarzWektor(b, N);
	usunMacierz(matrix, N);
	cout << endl;

	return 0;
}


double **stworzMacierz(int n) {
	double **matrix;
	matrix = new double *[n];
	for (int i = 0; i < n; i++)
		matrix[i] = new double[n];
	return matrix;
}

void uzupelnijMacierz(double **matrix, int n) {
	double tmp[4][4] = { { 1.0, 20.0, -30.0, -4.0 },
	{ 4.0, 20.0,-6.0, 50.0 },
	{ 9.0, -18.0, 12.0, -11.0 },
	{ 16.0, -15.0, 14.0, 130.0 } };
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			matrix[i][j] = tmp[i][j];
}

void usunMacierz(double **matrix, int n) {
	for (int i = n - 1; i >= 0; i--)
		delete[]matrix[i];
	delete[]matrix;
}
void pokarzWektor(double *b, int n)
{
	for (int i = 0; i < n;i++)
		cout << b[i] << endl;
}
void pokarzMacierz(double **matrix, int n) {
	//cout.setf(ios::basefield);
	//cout.precision(3);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
			cout <<matrix[i][j] << "\t";
		cout << endl;
	}
}

// częściowy wybór elementu podstawowego w przypadku, gdy dana wartość na przekątnej macierzy jest równa 0
int wybor_elementu_podstawowego(double **matrix, int n, int j, int *index) {
	int row;
	for (int i = n; i < j; i++) {
		if (fabs(matrix[index[i]][n]) < fabs(matrix[index[i + 1]][n]))
			row = index[i + 1];
		else
			row = index[i];
	}
	return row;
}

void gauss(double **matrix, int *index) {
	int row;
	double v;
	for (int k = 0; k < 3; k++) {
		if (matrix[index[k]][index[k]] == 0.0) {
			row = wybor_elementu_podstawowego(matrix, index[k], 3, index);
			index[row] = index[k];                                                     // zapisuje w tablicy i podmieniam indeksy
			index[k] = row;
		}
		for (int i = k + 1; i < 4; i++) {
			v = matrix[index[i]][k];
			for (int j = k + 1; j < 4; j++)
				matrix[index[i]][j] = matrix[index[i]][j] - matrix[index[k]][j] * (v / matrix[index[k]][k]);	//wyliczanie wartości macierzy U z eliminacji Gaussa
			matrix[index[i]][k] = v / (matrix[index[k]][k]);                     // Zapisanie L w tej samej macierzy co U
		}																			// L macierz współczynników przez które mnożymy kolejne wiersze (a[i][k]/a[k][k])
	}
}

// Ly=b  , wyznaczenie y
void macierzL(double **l, double *b, int *index, int n) {
	double sum = 0;
	for (int i = 0; i <= n; i++) {
		for (int j = 0; j < i; j++)
			sum = sum + l[index[i]][j] * b[index[j]];
		b[index[i]] = (b[index[i]] - sum) / 1.0;			//macierz L posiada jedynki na głownej przekątnej
		sum = 0;
	}
}

// Ux=y , wyznaczenie x
void macierzU(double **u, double *b, int *index, int n) {
	double sum = 0;
	for (int i = n; i >= 0; i--) {
		for (int j = i + 1; j <= n; j++)
			sum = sum + u[index[i]][j] * b[index[j]];
		b[index[i]] = (b[index[i]] - sum) / (u[index[i]][i]);
		sum = 0;
	}
}