function [L,U,p] = lutx(A)

%*****************************************************************************80
%
%% lutx() is a textbook version of the triangular factorization algorithm.
%
%  Discussion:
%
%    [L,U,p] = lutx(A) produces a unit lower triangular matrix L,
%    an upper triangular matrix U, and a permutation vector p,
%    so that L*U = A(p,:)
%
%  Licensing:
%
%    Copyright 2014 The MathWorks, Inc.
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625,
%    ebook: https://www.mathworks.com/moler/chapters.html
%
[n,n] = size(A);
p = (1:n)';

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   [r,m] = max(abs(A(k:n,k)));
   m = m+k-1;

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)
   
      % Swap pivot row
      if (m ~= k)
         A([k m],:) = A([m k],:);
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      i = k+1:n;
      A(i,k) = A(i,k)/A(k,k);

      % Update the remainder of the matrix
      j = k+1:n;
      A(i,j) = A(i,j) - A(i,k)*A(k,j); 
   end
end

% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);
