  function [L,d,M] = LDMinc(A,p,q)
% A is nxn and has an LDM factorization
% L and M are nxn and unit lower triangular and d is a column n vector 
% such that A = LDM' where D = diag(d).
[n,n] = size(A);
d = zeros(n,1);
% Pre-allocate storage for the triangles..
NL = nnz(tril(A))+n*p;
L = sparse([],[],[],n,n,NL)+speye(n,n);
NM = nnz(triu(A))+n*q;
M = sparse([],[],[],n,n,NM)+speye(n,n);
for j=1:n
    % Compute d(k), the kth column of L and the kth column of M...
    s    = d(1:j-1).*M(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;                 
    end
    L(j+1:n,j) = v(j+1:n)/d(j);
    
    t    = d(1:j-1).*L(j,1:j-1)';
    w(j+1:n) = A(j,j+1:n) - t'*M(j+1:n,1:j-1)';
    
    % Allocation for M(j+1:n,j)...
    m = nnz(A(j,j+1:n))+q;
    if m<n-j
        % We have to force some entries in w(j+1:n) to zero so determine 
        % where the unimportant nonzeros in w(j+1:n) are...
        inz     = find(w(j+1:n))+j;                 
        % inz = vector of nonzero component indices
        [y,idx] = sort(abs(w(inz)),'descend');    
        % inz(idx(m+1:end)) = vector of unimportant component indices
        w(inz(idx(m+1:end))) = 0;                 
    end
    M(j+1:n,j) = w(j+1:n)/d(j);
end

% same deduction ideas as in LDLinc 