% Naive Gaussian elimination
%
% Input:   Matrix A and vector b for linear system Ax = b.
%
% Outputs: Upper triangular matrix U and vector c for equivalent
%          system Ux = c.
%
% Note:    The computed matrix is only "virtually" upper triangular,
%          since 0's are not actually written in those positions
%          below the main diagonal.
%
function [U, c] = naiveGE(A, b)  
    % Ensure that A is square
    [n, m] = size(A);
    if n ~= m
        error('Square matrix expected.');
    end
    
    % Create a new column vector
    c = b;
    
    % Make a copy of A
    U = A;
    
    % Perform naive elimination...
    for j = 1 : n - 1
        % Slide down the main diagonal, introducing 0's below the diagonal
        if abs(U(j, j)) < eps
            error('Zero or near-zero pivot encountered.')
        end
        
        % Process each of the rows below the element A(j, j)
        for i = j + 1: n
            % Determine the multiplier needed to introduce a 0 in row i
            mult = U(i, j)/U(j, j);
            
            % Output multiplier information
            fprintf('Row(%d) = Row(%d) - (%8.5f) * Row(%d)\n', i, i, mult, j);
            
            % Perform the row operation needed to introduce the zero
            % ... Row i = Row i - mult * Row j
            for k = j + 1 : n
                U(i, k) = U(i, k) - mult * U(j, k);
            end
            
            % Adjust the right-hand side using the same multiplier
            c(i) = c(i) - mult * c(j);
        end
    end
end