function [x, U, bhat] = megp0 (A, b)
%
% function [x, U, bhat] = megp0 (A, b)
%
% Eliminazione di Gauss con pivoting (scambio delle righe).
n = length (A);
A = [A, b];
for k = 1:n - 1
  [~, idx] = max (abs (A(k:n, k)));
  temp = A(k, :);
  A(k, :) = A(idx + k - 1, :);
  A(idx + k - 1, :) = temp;
  for i = k + 1:n
    mik = A(i, k) / A(k, k);
    A(i, k) = 0;
    for j = k + 1:n + 1
      A(i, j) = A(i, j) - mik * A(k, j);
    end
  end
end
U = A(:, 1:n);
bhat = A(:, n + 1);
x = bs (U, bhat);
%!test
%! A = randn (4);
%! b = randn (4, 1);
%! x = megp0 (A, b);
%! assert (x, A\b, 1e-12)
%!test
%! A = [0,0,1;0,1,1;1,1,1];
%! b = [1;1;1];
%! x = megp0 (A, b);
%! assert (x, A\b, 1e-12)
