  
function P1Grade()
% Checks out EvecProd
% Generate a problem...
n = 1000;
S = triu(rand(n,n))+n*diag(rand(n,1));
T = triu(rand(n,n))+n*diag(rand(n,1));
k = 999;
clc

tic
xYOU = EvecProd(S,T,k);
tYOU = toc;
errYOU = norm(S*(T*xYOU) - S(k,k)*T(k,k)*xYOU);

tic
xCVL = EvecProdCVL(S,T,k);
tCVL = toc;
errCVL = norm(S*(T*xCVL) - S(k,k)*T(k,k)*xCVL);

fprintf('|| x ||                     = %10.6e\n',norm(xYOU))
fprintf('||S*T*x - S(k,k)*T(k,k)*x|| = %10.6e\n',norm(S*T*xYOU - S(k,k)*T(k,k)*xYOU))
fprintf('tYOU/tCVL = %10.3f\n',tYOU/tCVL)





function x = EvecProdCVL(S,T,k)
  % S and T are nxn upper triangular matrices and k is an integer with 1<=k<=n.
  % Assume that the eigenvalues of ST are distinct, i.e., S(1,1)T(1,1),...,S(n,n)T(n,n)
  % are distinct.
  % x is a unit 2-norm column vector so that S*T*x = lambda*x  where lambda = S(k,k)T(k,)

[n,n] = size(S);
x = zeros(n,1);
x(k) = -1;
if k==1
    return
else
    % Solve (U11-lambda*I)x(1:k-1) = b(1:k-1) where
    %     U11 = S(1:k-1,1:k-1)*T(1:k-1,1:k-1), 
    b(1:k) = S(1:k,1:k)*T(1:k,k);
    lambda = S(k,k)*T(k,k);
    x(k-1) = b(k-1)/(S(k-1,k-1)*T(k-1,k-1)-lambda);
    wc = zeros(k-1,1);
    wc(k-1) = T(k-1,k-1)*x(k-1);
    for j=k-2:-1:1
        % Compute x(j) using P3.1.5 ideas...
        g = T(j,j+1:k-1)*x(j+1:k-1);
        x(j) = (b(j) -S(j,j)*g - S(j,j+1:k-1)*wc(j+1:k-1))/...
                  (S(j,j)*T(j,j)-lambda);
        wc(j) = g+T(j,j)*x(j);
    end
    x = x/norm(x);
end
    
