  function [nNeg,nPos] = Sylvester(A)
% A is an nxn nonsingular symmetric matrix
% nNeg and nPos are the number of A's eigenvalues that are
% negative and positive.
[n,n] = size(A);
% P'*A*P = LDL'
[L,D,p] = ldl(A,'vector');
% Process the 2x2 blocks...
k=1;
while k<=n-1
    if abs(D(k+1,k))>0
        [Qk,Dk] = schur(D(k:k+1,k:k+1));
        D(k:k+1,k:k+1) = Dk;
        k = k+2;
    else
        k = k+1;
    end
end
nPos = sum(diag(D)>0);
nNeg = sum(diag(D)<0);

% -1 for using find excessively
% -3 for applying schur to D