% Eliminacja Gaussa z częściowym wyborem elementu głównego
clc; clear

% Oznaczony___________________________________
%a = [2 1 4 -2; -3 4 2 -1; 3 5 -2 1; -2 3 2 4];
%b= [19; 1; 8; 13];

a = [-3 1 1 -1; 3 6 -2 1; 2 0 5 -2; -2 3 2 9];
b = [2; -19; 25; -38];



% Nieoznaczony__________________________________

%a = [2 1 4 -2; -3 4 2 -1; -4 -2 -8 4; -2 3 2 4];
%b= [19; 1; -38 ; 13];


% Sprzeczny__________________________________

%a = [2 1 4 -2; -3 4 2 -1; -4 -2 -8 4; -2 3 2 4];
%b= [19; 1; -30 ; 13];


% Metoda _________________________________________________________________

    
    disp('Macierz główna.')
    disp('----------------------------')
    disp(a)

    [n, m] = size(a);

    if n ~= m
        error('Błąd: "a" musi być macierzą kwadratową.')
    end

    a = [a, b]; % Tworzymy macierz uzupełnioną

    
    disp('Macierz uzupełniona.')
    disp('----------------------------')
    disp(a)

    % Zaczyna się proces eliminacji

    disp('Eliminacja w przód.')
    disp('----------------------------')
    for i = 1:(n - 1)

        % Wybór wiersza głównego
        for j = (i + 1):n

            if abs(a(j, i)) > abs(a(i, i))
                a([i j], :) = a([j i], :); % Zamiana wierszy
            end

        end

        for j = (i + 1):n
            a(j, :) = a(j, :) - a(i, :) * (a(j, i) / a(i, i));
  
            % Sprawdzenie czy mamy wiersz zerowy
            licz = 0;
            for p = 1:m+1
                if a(j, p) == 0
                    licz = licz + 1;
                end
            end
            
            if licz == m+1
                error('Układ jest nieoznaczony');
            end

            if (licz == m && a(n,m+1) ~= 0)
                error('Układ jest sprzeczny');
            end
        end

            disp(a);
    end

    disp('Eliminacja w tył.')
    disp('----------------------------')

    for i = n:-1:2
        for j = i-1:-1:1
            a(j, :) = a(j, :) - a(i, :) * (a(j, i) / a(i, i));
        end
        disp(a);
    end

    % Generowanie rozwiązania

    for i = 1:n
        x(i) = a(i,n+1)/a(i,i);
    end

    disp('Rozwiązanie.')
    disp('----------------------------')
    disp('X = ')
    disp(x)