function [L, U, piv] = GEpiv(A);
%
% Pre:
% A n-by-n
%
% Post:
% L n-by-n unit lower triangular with |L(i, j)|<=1.
% U n-by-n upper triangular
% piv integer n-vector that is a permutation of 1:n.
%
% A(piv, :) = LU
[n, n] = size(A);
piv = 1:n;
for k=1:n-1
[maxv, r] = max(abs(A(k:n, k)));
q = r+k-1;
piv([k q]) = piv([q k]);
A([k q], :) = A([q k], :);
if A(k, k) ~= 0
A(k+1:n, k) = A(k+1:n, k)/A(k, k);
A(k+1:n, k+1:n) = A(k+1:n, k+1:n) - A(k+1:n, k)*A(k, k+1:n);
end
end
L = eye(n, n) + tril(A, -1);
U = triu(A);