function [x,l] = gauss_piv(A,b)
% Gauss Elimination with partial pivoting. 
% Output: The solution x
%         The pivot vector l
%
n = length(b);
x = zeros(n,1);
l = [1:n]';

% Forward elimination
for k=1:n-1
    [tmp,ind] = max(abs(A(l(k:end),k)));
    r = k-1+ind;
    ncopy = l(k);
    l(k)=l(r);
    l(r)= ncopy; 
    
    for i=k+1:n
        m = A(l(i),k)/A(l(k),k);
        A(l(i),k)=0;
        for j=k+1:n
            A(l(i),j) = A(l(i),j)-m*A(l(k),j);
        end
        b(l(i)) = b(l(i))-m*b(l(k));
    end
end

% backward solution
x(n) = b(l(n))/A(l(n),n);
for k=n-1:-1:1
    res = b(l(k));
    for s=k+1:n
        res = res-A(l(k),s)*x(s);
    end
    x(k) = res/A(l(k),k); 
end


