function [L,U,p] = myLU_pivot(A)
% Call : function [L,U,p] = myLU_pivot(A)
% Input : A - nxn matrix
% Output: L - nxn unit lower triangular
%         U - nxn upper triangular
%         p - vector of length (n-1),
%             p(j) = k means at step j rows j&k are swapped


[~,n] = size(A);

L = eye(n);
for i = 1:n-1
    
    % find the largest entry in |.| on the ith column below (i,i)
    [~,k] = max(abs(A(i:n,i)));
    p(i) = i + (k-1);
    
    % swap rows i & p(i)
    temp = A(p(i),i:n);
    A(p(i),i:n) = A(i,i:n);
    A(i,i:n) = temp;
    
    % swap also the multipliers at rows i & p(i)
    temp = L(p(i),1:i-1);
    L(p(i),1:i-1) = L(i,1:i-1);
    L(i,1:i-1) = temp;
    
    for j = i+1:n
        L(j,i) = A(j,i)/A(i,i);
        A(j,i:n) = A(j,i:n) - L(j,i)*A(i,i:n);
        A(j,i) = 0;
    end
end

U = A;
return;

     

