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)
%
%% initialization and preparation
% check if the dimensions of the problem somehow match and return an error
% message if not
dimA = size(A);
if(length(dimA) > 2 || dimA(1)~= dimA(2))
    error('A is not a aquare matrix.')
end

n = dimA(2);

% initialize the index vector l
l = 1:n;

% initialize the vector of ratios needed for determining the pivot element
r = zeros(n,1);

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

%% forward elimination
for k=1:n-1
    % compute the ratios between (relevant) entries in column scale and
    % respective scale.
    r(k:end) = A(l(k:end),k)./s(l(k:end));
    % store the position where the maximum is attained in the variable pos.
    % Note that the position is calculated in the truncated vector
    % containing only the entries from k to n. Thus the actual position is
    % obtained by adding k-1.
    [temp pos] = max(r(k:end));
    % If necessary simulate a row permutation.
    if pos > 1
        l([k pos+k-1]) = l([pos+k-1 k]);
    end
    % compute the multipliers and store them in the matrix A
    A(l(k+1:end),k) = A(l(k+1:end),k)/A(l(k),k);
    
    % subtract multiples of line k from the next lines
    A(l(k+1:end),k+1:end) = A(l(k+1:end),k+1:end)-A(l(k+1:end),k)*A(l(k),k+1:end);
end