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

double EPS = 1.0e-10;
const int SIZE = 4;

void printMatrix(double arr[][SIZE]) {
	for(int i=0;i<SIZE;i++)
	{
		for (int j = 0; j < SIZE; j++)
		{
			std::cout << arr[i][j] << "  ";
		}
		std::cout << std::endl;
	}
	std::cout << std::endl;
}
void printVector(double arr[]){
	for (int i = 0; i < SIZE; i++) {
		std::cout << arr[i] << std::endl;
	}
	std::cout << "\n\n";
}
void printVectorB(double* b) {
	for (int i = 0; i < SIZE; i++) {
		std::cout << *(b + i) << std::endl;
	}
	std::cout << "\n\n";
}
int max(double arr[][SIZE], int kol){
	double max = 0;
	int maxID = 0;
	for (int i = 0; i < SIZE; i++)
	{
		if (fabs(arr[i][kol]) > max)
		{
			max = fabs(arr[i][kol]);
			maxID = i;
		}
	}
	return maxID;
}
double* opV(double wektor[], double piv[][SIZE], int swaps[][2]){
	static double vect[SIZE];
	for (int i = 0; i < SIZE; i++)
	{
		vect[i] = wektor[i];
	}
	for (int i = 0; i < SIZE; i++)
	{
		while (swaps[i][0] != swaps[i][1]) {
			std::swap(vect[swaps[i][0]], vect[swaps[i][1]]);
			std::cout << "Zamieniam " << swaps[i][0] << " z " << swaps[i][1] << std::endl;
			break;
		}
	}
	//std::cout << "Po zamianie wierszy: ";
	//printVector(vect);
	for (int k = 0; k < SIZE; k++)
	{
		for (int i = 1; i < SIZE; i++)
		{
			if (k < i){
				vect[i] -= (vect[k]*piv[i][k]);
			}
		}
	}
	return vect;
}
int main() {
	//Napisz program w języku „C/C++”, realizujący dekompozycję LU macierzy A, przy zastosowaniu
	//eliminacji Gaussa z częściowym wyborem elementu podstawowego, a następnie rozwiązujący układ
	//równań Ax = b.
																			//
	double A[SIZE][SIZE] = {{1.,20.,-30.,-4.}, {4.,20.,-6.,50.}, {9.,-18.,12.,-11.}, {16.,-15.,14.,130.}};
	double B[SIZE] = {0. , 114., -5., 177.};

	//double A[SIZE][SIZE] = { {6.,-2.,2.,4.}, {12.,-8.,6.,10.}, {3.,-13.,9.,3.}, {-6.,4.,1.,-18.} };
	//double B[SIZE] = { 12.,34.,27.,-38. };

	double L[SIZE][SIZE] = { {1.,0.,0.,0.},{0.,1.,0.,0.},{0.,0.,1.,0.},{0.,0.,0.,1.} };
	double X[SIZE] = {};
	int zamiany[SIZE][2] = {{0,0},{0,0},{0,0},{0,0}};
	int id = 0;
	printMatrix(A);
	//A = L * U
	//1.operacja na macierzy A
	for (int k = 0; k < SIZE; k++)//kolumny
	{
		for (int i = 1; i < SIZE; i++)//wiersze
		{
			if (k < i) // elementy pod przekątną
			{
				if (A[k][k] == 0.) // jezeli el na przekatnej == 0 to nie mozemy przez niego podzielic
				{
					int maxIndex = max(A, k);
					double wiersz[SIZE];

					for (int i = 0; i < SIZE; i++)
					{
						wiersz[i] = A[maxIndex][i];	//zapisanie wartosci wiersza do zamiany
					}

					for (int i = 0; i < SIZE; i++)	//zamiana wiersza
					{
						A[maxIndex][i] = A[k][i];
						A[k][i] = wiersz[i];
					}
					zamiany[id][0] = k;
					zamiany[id][1] = maxIndex;
					std::cout<< "Do zamiany wyrazy o indexach " << zamiany[id][0] << " z " << zamiany[id][1]<<"\n" << std::endl;
					id++;
				}
				double pivot = A[i][k] / A[k][k];
				L[i][k] = pivot;	//mnoznik
				for (int col = 0; col < SIZE; col++)//kolumny
				{
				A[i][col] -= (A[k][col] * pivot);
				}
				//B[i] -= (B[k] * pivot);
			}
		}
	}
	printMatrix(A);
	printMatrix(L);
	//operacja na wektorze b
	double* b = opV(B, L, zamiany);
	//std::cout << "B po operacjach wierszy: ";
	printVectorB(b);
	//wyliczenie x
	for (int i = SIZE - 1; i >= 0; --i) {
		double temp = 0;
		for (int j = i + 1; j < SIZE; ++j) {
			temp += A[i][j] * X[j];
		}
		X[i] = ((*(b+i)) - temp) / A[i][i];
	}
	printVector(X);
	getchar();
}