#pragma warning (disable:4996)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define R 7 

main()
{
	double A[R][R] = { 0 }, X[R] = { 0 }, B[R] = { 0 }, L[R][R] = { 0 }, D[R][R] = { 0 }, U[R][R] = { 0 }, Z[R] = { 0 }, Y[R] = { 0 };//macierze deklaracja
	int i, j, k, rozmiar,pozycja,n;//inty    n jest tylko do wektora x wykorzystany
	double iloczyn = 1, suma = 0, ysuma = 0, yiloczyn = 1;//do obliczania potrzebne
														  /*Otwieranie pliku*/
	FILE *plik1;
	plik1 = fopen("plik01.txt", "r");
	if (plik1 == NULL)
	{
		printf("PLIK 1 JEST PUSTY LUB NIE ISTNIEJE\n");
		system("pause");
	}
	else
	{
		/*KOMUNIKAT ZE PLIK ISTNIEJE*/
		printf("PLIK 1 ISTNIEJE :) \n");
		fscanf(plik1, "%d", &rozmiar);					//WCZYTANIE ROZMIARU MACIERZY
		printf("Rozmiar macierzy: %d\n\n", rozmiar);
		for (i = 0; i < rozmiar; i++)					//wczytujemy macierz A 
		{
			for (j = 0; j < rozmiar; j++)
			{
				fscanf(plik1, "%lf", &A[i][j]);
			}
		}
		//wczytanie wektora B
		for (i = 0; i < rozmiar; i++)
		{
			fscanf(plik1, "%lf", &B[i]);
		}


		/*ROZKLAD LDU ROZKLAD LDU ROZKLAD DLU*/
		/*KROK 1. DIAGONALKA I WIERSZ/KOLUMNA*/
		D[0][0] = A[0][0];
		for (i = 0; i < rozmiar; i++)					//L na kolumnie U na wierszu
		{
			L[i][0] = A[i][0] / D[0][0];
			U[0][i] = A[0][i] / D[0][0];
		}
		//TO JEST OK WIEC NIE RUSZAJ	
		/**/
		for (i = 1; i <rozmiar; i++)
		{
			L[i][i] = 1;
			U[i][i] = 1;
		}
		/*KROK 2. RESZTA DO OBLICZANIA MASAKRA*/
		for (i = 1; i < rozmiar; i++)
		{
			for (j = 1; j < rozmiar; j++)
			{
				//DIAGONALNA gdy i=j
				if (i == j)
				{
					suma = 0; iloczyn = 1;
					for (k = 0; k < i; k++)
					{
						iloczyn = L[i][k] * U[k][i] * D[k][k];
						suma += iloczyn;
					}
					D[i][i] = A[i][i] - suma;		//wynik diagonalki
				}
				//maczierz L gdy i>j
				else if (i > j)
				{
					suma = 0; iloczyn = 1;
					for (k = 0; k < j; k++)
					{
						iloczyn = L[i][k] * U[k][j] * D[k][k];
						suma += iloczyn;
					}
					L[i][j] = (A[i][j] - suma) / D[j][j];
					/*MACIERZ U TERAZ*/
					k = i;//zamiana wspolrzednych
					i = j;
					j = k;
					suma = 0; iloczyn = 1;
					for (k = 0; k < j; k++)
					{
						iloczyn = L[i][k] * U[k][j] * D[k][k];
						suma += iloczyn;
					}
					U[i][j] = (A[i][j] - suma) / D[i][i];
					k = j;
					j = i;
					i = k;
				}
			}
		}
		/* wypisanie LDU */
		printf("Macierz L \n");
		for (i = 0; i < rozmiar; i++)
		{
			for (j = 0; j < rozmiar; j++)
			{
				printf("%lf\t", L[i][j]);
			}
			printf("\n");
		}
		printf("\nMacierz D  \n");
		for (i = 0; i < rozmiar; i++)
		{
			for (j = 0; j < rozmiar; j++)
			{
				printf("%lf\t", D[i][j]);
			}
			printf("\n");
		}
		printf("\nMacierz U \n");
		for (i = 0; i < rozmiar; i++)
		{
			for (j = 0; j < rozmiar; j++)
			{
				printf("%lf\t", U[i][j]);
			}
			printf("\n");
		}
		/*WYPISANIE MACIERZY A I WEKTORA X*/
		/*
		printf("\nMacierz A \n");
		for ( i = 0; i < rozmiar; i++)
		{
		for (j = 0; j < rozmiar; j++)
		{
		printf("%lf  ", A[i][j]);
		}
		printf("\n");
		}
		printf("\nWektor X \n");
		for (i = 0; i < rozmiar; i++)
		{
		printf("%lf  ", X[i]);
		}
		*/
		
		/*Wektor Y */
		Y[0] = B[0];
		ysuma = 0, yiloczyn = 1;
		for ( i = 1; i < rozmiar; i++)
		{
			for ( j = 0; j < rozmiar; j++)
			{
				yiloczyn = L[i][j] * Y[j];
				ysuma += yiloczyn;
			}
			Y[i] = B[i] - ysuma;
			ysuma = 0, yiloczyn = 1;
		}
		/*Wektor Z */
		for ( i = 0; i < rozmiar; i++)
		{
			Z[i] = Y[i] / D[i][i];
		}
		/*Wektor X */
		pozycja = rozmiar - 1;
		X[pozycja] = Z[pozycja];
		for ( i = pozycja; 0 < i ; i--)
		{
			n = i - 1;
			for ( j = i; j < rozmiar; j++)
			{
				yiloczyn = U[n][j] * X[j];
				ysuma += yiloczyn;
			}
			X[n] = Z[n] - ysuma;
			ysuma = 0, yiloczyn = 1;
		}

		printf("\n\nWEKTOR XXXXXXX\n");
		for ( i = 0; i < rozmiar; i++)
		{
			printf("%lf   ", X[i]);
		}

		printf("\n\n");
		system("pause");
	}
}