function [x] = gem(A, b, piv)
%GEM  Gaussian elimination
%   X=GEM(A,B,PIV) solves the square linear system
%   AX=B using Gaussian Elimination without pivoting
%   if PIV = 0 or with partial pivoting if PIV = 1.
%
%   X = GEM(A,B) as before assuming PIV = 1.
%

% set piv to 1 if it is not specified as an input
if nargin<3
    piv = 1;
end

% find the size of the linear system
[m, n] = size(A);

% begin the solution process
for k=1:n-1
    % do k-th step of the elimination
    if piv
        % search for the element of maximun modulus
        % along k-th column and, if it find in position
        % A(i,k) with i>k, exchange rows k and i 
        [ maxPiv, iPiv ] = max( abs( A(k:n,k) ) );
        % index iPiv is referred to the length of the
        % vector A(k:n,k) so we have to find the proper
        % value of i in the matrix A
        i = k+iPiv-1;
        
        % exchange k-th row with i-th row if necessary
        if maxPiv>abs( A(k,k) )
            % exchange A elements
            savekthRow = A(k,k:n);
            A(k,k:n) = A(i,k:n);
            A(i,k:n) = savekthRow;
            % exchange b elements
            savekthb = b(k);
            b(k) = b(i);
            b(i) = savekthb;
        end
    end
        
    % update the right bottom part of the matrix A
    for i=k+1:n
        % update i-th row of A
        lik = A(i,k)/A(k,k); % multiplier
        A(i,k+1:n) = A(i,k+1:n)-lik*A(k,k+1:n);
        % update i-th element of b
        b(i) = b(i) - lik*b(k);        
    end
    
    % set to zeros the element in the k-th column
    % below A(k,k)
    A(k+1:n,k) = 0;
end

% backward substitution method
x = zeros(n,1);
x(n) = b(n)/A(n,n);
for k=n-1:-1:1
    % compute x(k)
    x(k) = ( b(k)-A(k,k+1:n)*x(k+1:n) )/A(k,k);
end
