%-------------------------------------------------------------------------------
%
%  This function computes a LU factorisation of matrix A. Given A, it 
%  outputs matrices L (lower triangular) and U (upper triangular)
%
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------

function [L,U]  =  MyLUFactorisation(A)
%--------------------------------  
% Check dimensions of A
%--------------------------------  
n  =  size(A,1); % número de filas de A 
m  =  size(A,2); % número de columnas de A
if abs(n-m)>0
  error('Chequee las dimensiones de A. La matriz debe ser cuadrada')  
end
%------------------------------------
% Inicializamos las matrices L y U con ceros. Trabajamos en formato full
% Implementación de la factorización en format sparse es difícil
%------------------------------------
L                      =  zeros(n);
U                      =  zeros(n);
for j=1:n  %  bucle sobre el número de columnas
    %  Introducimos la diagonal de L: L(j,j)=1
    L(j,j)             =  1;
    for i=1:n  %  bucle sobre el número de filas
        %-----------------------------------
        % Situation where i>j
        %-----------------------------------
        if i>j
           L(i,j)      =  A(i,j);
           for k=1:j-1
               L(i,j)  =  L(i,j) - L(i,k)*U(k,j);    
           end 
           L(i,j)      =  L(i,j)/U(j,j);
        end
        %-----------------------------------
        % Situation where i<=j
        %-----------------------------------
        if i<=j
           U(i,j)      =  A(i,j);
           for k=1:i-1
               U(i,j)  =  U(i,j) - L(i,k)*U(k,j);    
           end 
        end
    end
end
  
end
