sdchol(A)

function T = sdchol(A)
%SDCHOL Cholesky-like decomposition for positive semidefinite matrices
%   T = SDCHOL(A) uses only the diagonal and upper triangle of A.  If
%   A is positive semi-definite, then T=SDCHOL(A) produces a matrix T
%   such that A = T'*T.  If A is positive definite, then T is the
%   upper triangular Cholesky factor returned by chol(A).  If A is not
%   positive definite, T may be neither square nor triangular.  If
%   A has negative eigenvalues, an error message will be printed.

    
[T,p] = chol(A);

if p > 0
    [m,n] = size(A);
    assert(m == n, 'Matrix must be square.');

    % It's too bad matlab's ldl() doesn't have an option forcing the
    % diagonal matrix to be diagonal (and possibly complex) instead
    % of block diagonal, as it'd probably be the fastest and most
    % stable algorithm here. But we'll make do with the Schur
    % decomposition instead.
    A = triu(A) + triu(A,1)'; % use only upper triangle
    [U, D] = schur(A, 'complex'); % A = U*D*U' with U unitary, D
                                  % diagonal (since A symmetric)
    D = diag(D);
    
    if ~isreal(D) 
        if any(imag(D) > 10*eps)
            error('Matrix must be positive semi-definite');
        else
            D = real(D);
        end
    end

    if any(D < 0)
        if any(D < 10 * eps)
            error('Matrix must be positive semi-definite');
        else
            D(D < 0) = 0;
        end
    end
    
    T = diag(sqrt(D)) * U';
end