function [L,U,P]= LU_factorization(A)
% This function returns the LU factorization of a matrix A
% By using Gauss elimination method with pivoting by row
% 
%       P*A = L*U     
% 
% Input  -> A (n x n matrix)
% Output -> L lower triangular matrix
%           U Upper triangular matrix
%           P is the permutation matrix
%
%     

N=size(A,1) ;  % size of the square matrix A
L=eye(N);      % here we initialize the matrix L to be the identity matrix
P=eye(N);      % here we initialize the permutation matrix

I=linspace(1,N,N)'; % with this vector we keep track of all permutations
                    % we do during Gauss elimintation (pivoting by row)

for i=1:N-1 

% Pivoting step
    [aaa,idx]=max(abs(A(i:end,i)));

    idx=idx+(i-1);
    if idx~=i
        b       = A(idx,:); % This is the idx-row of the matrix A
        A(idx,:)= A(i,:);
        A(i,:)  = b;
        
        a   = I(i);
        I(i)=I(idx);
        I(idx)=a;
    end

% Elimination step
    for j=i+1:N
        m=A(j,i)/A(i,i);
        A(j,i:end)=A(j,i:end)-m*A(i,i:end);        
        A(j,i)=m;                         
    end
    A
    I
end

L=eye(N)+tril(A,-1);
U=triu(A);
P=P(I,:);

