%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
%
%  This function computes a Cholesky factorisation of matrix A. Given A, it 
%  outputs matrices L, such that A =  L*L'
%
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------

function L     =  MyCholeskyFactorisation(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 la matriz
% con ceros. Trabajamos en formato full
% Implementación de la factorización 
% en format sparse es muy difícil
%------------------------------------
L                      =  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)             =  A(j,j); 
    for k=1:j-1
        L(j,j)         =  L(j,j) - L(j,k)^2;
    end
    L(j,j)             =  sqrt(L(j,j));     
    for i=1:n  %  bucle sobre el número de filas
        if abs(i-j)>0
           L(i,j)           =  A(i,j);
           for k=1:j-1
               L(i,j)         =  L(i,j) - L(i,k)*L(j,k);
           end
           L(i,j)           =  L(i,j)/L(j,j);
        end
    end     
end 
  
end