function x = gauss_elim_vector_piv(A,b)

% A should be square
% size(A) = [length(x) length(x)]

n = size(A,1);

p = 1:n;
for k = 1:(n-1)
    rows = k:n;
    [y,v] = max(abs(A(p(rows),k)));
    i = rows(v);
    p([k i]) = p([i k]);
    for i = (k+1):n
        mult = A(p(i),k) / A(p(k),k);
        A(p(i),k:n) = A(p(i),k:n) - mult * A(p(k),k:n);
        b(p(i)) = b(p(i)) - mult * b(p(k));
    end
end

x = zeros(n,1);
for k = n:-1:1
    sum = b(p(k)) - A(p(k),(k+1):n) * x((k+1):n);
    x(k) = sum / A(p(k),k);
end