function [Q,D,mu,x] = SchurPlus(A,v)
% A is an n-by-n symmetric positive definite matrix.
% v is a column n-vector
% Q'*A*Q = D is the Schur decomposition of A.
% mu is the largest eigenvalue of A + v*v' and x is an associated eigenvector, 
% normalized so  norm(x,2) = 1. Assume that the largest eigenvalue of A and
% A + v*v' are unique

n = length(v);
[Q,D] = schur(A);
% Note that Q'(A + vv')Q = D + ww'  where w = Q'v
w = Q'*v;
d = diag(D);
% Idea: Apply power method to (D + ww') to get dominant eigenvalue mu
% and associated eigenvector z: (D + ww')z = mu*z

% Initializations...
z = randn(n,1); z = z/norm(z);
[y,mu] = Prod(d,w,z);
tol = norm(A+v*v','inf')*10^-14;

% Note that if x is almost an eigenvector then y = (D+ww')x \approx mu*x
% where mu = x'*y.
while norm(y - mu*z)>tol
    z = y/norm(y);
    [y,mu] = Prod(d,w,z);
end

% If (D+ww')z = mu*z then (A+vv')x = mu x where
x = Q*z;

function [y,mu] = Prod(d,w,z)
% y = (diag(d) + ww')z
% mu = z'*y
y = d.*z + (w'*z)*w;
mu = z'*y;
