function P4Grade()
% Checks out SpecialEig
clc
n = 1000;
d = randn(n,1);
v = randn(n,1);
ev = diag(schur(diag(d)+v*v'));
ev = sort(ev);
[z,idx] = sort(abs(ev));
for k = [20 40 100 200 400]
   tic
   evSA = SpecialEig(d,v,k,'SA');
   tSA = toc;
   errSA = norm(evSA-sort(ev(1:k)));
   fprintf('n=%1d  k = %3d  errSA = %10.3e  tSA = %6.3f\n',n,k,errSA,tSA)
   tic
   evLA = SpecialEig(d,v,k,'LA');
   tLA = toc;
   errLA = norm(evLA-sort(ev(n-k+1:n)));
   fprintf('n=%1d  k = %3d  errLA = %10.3e  tLA = %6.3f\n',n,k,errLA,tLA)
   tic
   evSM = SpecialEig(d,v,k,'SM');
   tSM = toc;
   errSM = norm(evSM-sort(ev(idx(1:k))));
   fprintf('n=%1d  k = %3d  errSM = %10.3e  tSM = %6.3f\n',n,k,errSM,tSM)
   tic
   evLM = SpecialEig(d,v,k,'LM');
   tLM = toc;
   errLM = norm(evLM-sort(ev(idx(n-k+1:n))));
   fprintf('n=%1d  k = %3d  errLA = %10.3e  tLM = %6.3f\n\n',n,k,errLM,tLM)
end

function eValues = SpecialEigCVL(d,v,k,what)
% d is a column n-vector with distinct nonzero entries
% v is a column n-vector with nonzero entries
% k is a positive integer that satifies 1<=k<n
% Let A = diag(d) + vv'.
% If what equals 'SA', then eValues is a column k-vector comprised 
%         of the k algebraically smallest eigenvalues of A.
% If what equals 'LA', then eValues is a column k-vector comprised 
%         of the k algebraically largest eigenvalues of A.
% If what equals 'SM', then eValues is a column k-vector comprised 
%         of the k smallest eigenvalues of A in magnitude
% If what equals 'LM', then eValues is a column k-vector comprised 
%         of the k largest eigenvalues of A in magnitude.
% In each case, eValues(1) < eValues(2) < ... < eValues(k).

   opts.issym = 1;
   n = length(d);
   if strcmp(what,'SA')|| strcmp(what,'LA') || strcmp(what,'LM')
       eValues = eigs(@(x) MatVecProd(x,d,v),n,k,what,opts);
   else
       eValues = eigs(@(x) InvMatVecProd(x,d,v),n,k,what,opts);
   end
   eValues = sort(eValues);




function y = MatVecProd(x,d,v)
% y = (diag(d)+vv')x
y = d.*x + (v'*x)*v;

function y = InvMatVecProd(x,d,v)
% y = (diag(d)+vv')\x
u = v./d;
z = x./d;
y = z - ((v'*z)/(1+v'*u))*u;
