function [L, U, P] = lup0_kij (A)
%
% function [L, U] = lup0_kij (A)
% function [L, U, P] = lup0_kij (A)
%
% Fattorizzazione LU con pivoting (con scambio delle righe).
n = length (A);
p = (1:n);
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;
  temp = p(k);
  p(k) = p(idx + k - 1);
  p(idx + k - 1) = temp;
  for i = k+1:n
    A(i, k) = A(i, k) / A(k, k);
    for j = k+1:n
      A(i, j) = A(i, j) - A(i, k) * A(k, j);
    end
  end
end
U = triu (A);
I = eye (n);
L = tril (A) - diag (diag (A)) + I;
if (nargout == 2)
   L(p, :) = L;
else
  P = I(p, :);
end
%!test
%! A = randn (10);
%! [LL, UU] = lup0_kij (A);
%! [L, U] = lu (A);
%! assert (L, LL, 1e-14)
%! assert (U, UU, 1e-14)
%! [LL, UU, PP] = lup0_kij (A);
%! [L, U, P] = lu (A);
%! assert (P, PP, 1e-14)
