function [A,l] = Gauss(A)
% Gauss performs the forward elimination procedure on the coefficient 
% matrix A using Gaussian elimination with scaled partial pivoting.
%
% INPUT: A (Coefficient matrix, assumed invertible.)
%
% OUTPUT: A (The resulting matrix from the forward elimination procedure)
%         l (Index vector. If the rows of A are viewed in the order given
%            in l, the upper triangular part holds the transformed 
%            coefficient matrix and the strictly lower triangular part
%            holds the multipliers used in the elimination procedure)
%

% Find number of unknowns
n = length(A);

% Compute scale vector
s = max(abs(A),[],2);

% Initialize index vector
l = 1:n;

for k = 1:(n-1)
        
        % Find max scaled pivot
        [~,p] = max(abs(A(l(k:n),k))./s(l(k:n)));
        
        % Swap elements in the index vector if row swap is required
        if p ~= 1
            temp = l(k);
            l(k) = l(k+p-1);
            l(k+p-1) = temp;
        end
        
        % Compute and store multipliers in "lower" part of A
        A(l(k+1:n),k) = A(l(k+1:n),k)/A(l(k),k);
        
        % Perform the eliminations
        A(l(k+1:n),k+1:n) = A(l(k+1:n),k+1:n)-A(l(k+1:n),k)*A(l(k),k+1:n);
        
end

end