function [P,L,U] = permFac(A)
% permFac computes P, L and U for a given matrix A such that PA = LU.
%
% INPUT: A (Matrix to factorize, assumed invertible)
% 
%
% OUTPUT: P (Permutation matrix)
%         L (Unit lower triangular matrix)
%         U (Upper triangular matrix)
%

n = length(A);
s = max(abs(A),[],2);
% Perform Gaussian elimination with scaled partial pivoting.
l = 1:n;
for k = 1:(n-1)
        [~,p] = max(abs(A(l(k:n),k))./s(l(k:n)));
        if p ~= 1
            temp = l(k);
            l(k) = l(k+p-1);
            l(k+p-1) = temp;
        end
        mult = A(l(k+1:n),k)/A(l(k),k);
        A(l(k+1:n),k) = mult;
        A(l(k+1:n),k+1:n) = A(l(k+1:n),k+1:n)-mult*A(l(k),k+1:n);
end
% Create the permuation matrix directly from the index vector.
P = zeros(n);
P((1:n)+(l-1)*n) = ones(n,1);

% Create L and U from the result of the forward elimination procedure on PA
L = tril(A(l,:),-1);
L(1:n+1:n^2) = ones(n,1);

U = triu(A(l,:));



end