function x = FindNegVec(A)
% A is n-by-n, symmetric, and has exactly one negative eigenvalue
% x is a unit 2-norm vector so that x'*A*x < 0

[n,n] = size(A);
% Compute the factorization P'*A*P = LDL' and
% note that if z'*D*z < 0, L'y = z, and x = P*y, then
% x'*A*x = y'(P'*A*P)y = y'(L*D*L')y = z'*D*z < 0
[L,D,p] = ldl(A,'vector');
% So let's get z.
% Does D have a 2-by-2 block or not?
idx = find(abs(diag(D,1))>0);
if isempty(idx)
% D is diagonal)
k = find(diag(D)<0);
% D(k,k) < = 0
z = zeros(n,1); z(k) = 1;
else
% This is the sole 2-by-2 block...
k = idx+1;
Dk = D(k-1:k,k-1:k);
% Compute its eigenvectors and eigenvalues...
[V,Lambda] = eig(Dk);
if Lambda(1)<1
z = [zeros(k-2,1);V(:,1);zeros(n-k,1)];
else
z = [zeros(k-2,1);V(:,2);zeros(n-k,1)];
end
end
% Solve the upper triangular system L'w = z and note z(k+1:n) = 0.
% (OK to ignore since overall is O(n^3) and this is an O(n^2) shortcut.)
w = [L(1:k,1:k)'\z(1:k);zeros(n-k,1)];
x = w(p);
x = x/norm(x);