function output = CholSemidef(A,CholSemidef_Tol)
% Copyright (c) 2001-2009 by Ales Cerny
%-----------------------------------------------------------
% Cholesky decomposition for positive semidefinite matrices 
%-----------------------------------------------------------

%************************************************************************%
% chapter11sect9a.m - supplementary program to                            %
% Ales Cerny (2009) Mathematical Techniques in Finance (2nd ed.)         % 
% Princeton University Press http://press.princeton.edu/titles/9079.html %
%************************************************************************%
    
% This code is provided 'as-is', without any express or implied warranty.  
% 
% Permission is granted to anyone to use this code for any purpose, 
% subject to the following restrictions:
% 
% 1. The origin of this code must not be misrepresented; you must not
%    claim that you wrote the original code.
% 2. Modified code versions must be plainly marked as such, and must not 
%    be misrepresented as being the original code.
% 3. This notice may not be removed from any source distribution.

% NOTICE TO STUDENTS: To avoid accusations of plagiarism, if you use this 
% code or its modifications in assessed work you should prepend it with a 
% note stating: 
%   "This is the original/modified version of the code chapter11sect9a.m by 
%    Ales Cerny (2009), Mathematical Techniques in Finance (2nd ed.),  
%    Princeton University Press. The original version is available from 
%    http://www.martingales.info/mtfweb2".
% A similar acknowledgement should appear prominently inside your written 
% report.

[m,n] = size(A);

if (m ~= n)
    error('error CholSemidef: Input matrix not square');
elseif ( sum(sum((A-A')^2))>0)
    error('error CholSemidef: Input matrix not symmetric');
elseif (A(1,1)<= 0) 
    error('error CholSemidef: Variance of the first variable must be positive');
end
    
sig = zeros(n,n);
sig(1,1) = sqrt(A(1,1));
if n == 1
 output = sig;
else
    sig(2:n,1)=A(2:n,1)/sig(1,1);
    idx = 1;        % idx contains indices of all linearly independent random variables
    for i = 2 : 1 : n-1
        aux1= sig(i,idx); % select columns corresponding to linearly independent variables
        aux2=A(i,i)-aux1*aux1';
        if (aux2 < 0)
            disp(sprintf('problem with random variable no: %3.0f', i));
            error('error CholSemidef: Input matrix not a covariance matrix (not positive semidefinite)');
        elseif (aux2 > CholSemidef_Tol)  % this random variable is not redundant
            sig(i,i)=sqrt(aux2);
        end
            if i < n
                for k = i+1 : 1 : n   %(i+1,n,1);
                aux1 = sig(k,idx);
                aux2 = sig(i,idx);
                sig(k,i)=(A(k,i)- aux1*aux2')/sig(i,i);    
                end
            end
            idx = [idx i];    %add this variable to the list of linearly independent ones
    end
    output = sig(:,idx);  % return all rows, but only those columns corresponding to lin. indep. vars
end