  function [L,d] = LDLincCVL(A,p)
% A is an nxn and symmetric positive definite matrix
% L is unit lower triangular and d is a column n-vector such that
% A = LDL' where D = diag(d).
[n,n] = size(A);
d = zeros(n,1);
N = nnz(tril(A))+n*p;
L = sparse([],[],[],n,n,N)+speye(n,n);
for j=1:n
    % Compute d(k) and the k-th column of L...
    s          = d(1:j-1).*L(j,1:j-1)';
    d(j)       = A(j,j) - L(j,1:j-1)*s;
    v(j+1:n)   = A(j+1:n,j) - L(j+1:n,1:j-1)*s;
    % Allocation for L(j+1:n,j)...
    m = nnz(A(j+1:n,j))+p;
    if m<n-j
        % We have to force some entries in v(j+1:n) to zero so determine 
        % where the unimportant nonzeros in v(j+1:n) are...
        inz     = find(v(j+1:n))+j;                 
        % inz = vector of nonzero component indices
        [y,idx] = sort(abs(v(inz)),'descend');    
        % inz(idx(m+1:end)) = vector of unimportant component indices
        v(inz(idx(m+1:end))) = 0;                 
        % The following is slower and involves lengh(n-j) vectors...
        %   [y,idx] = sort(v(j+1:n),'descend');    
        %   v(idx(m+1:end)) = 0;
    end
    L(j+1:n,j) = v(j+1:n)/d(j);
end

% -3 or so  if your sort step involves O((n-j)) flops
% -4 or so if you work with full format L and then putting L = sparse(L) at the end 
