function [L,U,p] = modlutx(A)
%MODLUTX  Triangular factorization.  Modified lutx; this has bland for loops.
%   [L,U,p] = modlutx(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U, and a permutation vector p,
%   so that L*U = A(p,:)

[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)
         for j = 1:n
           temp = A(k,j);
           A(k,j) = A(m,j);
           A(m,j) = temp;
         end
         temp = p(k);
         p(k) = p(m);
         p(m) = temp;
      end

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

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

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