function P1Grade()

clc
m = 100;
n = 100;
[U,R] = qr(randn(m,n));
tol = 10^(-12);

% Example 1
alfa = -0.9;
beta = 0.0;
kYOU = nEigs(U,alfa,beta);
kCVL = nEigsCVL(U,alfa,beta);
ev = real(eig(U));
kLONG = sum(alfa-tol<=ev & ev<=beta+tol);
fprintf('\nm = %2d   n = %2d  alfa = %6.3f  beta = %6.3f\n',m,n,alfa,beta)
fprintf('kYOU  = %1d\nkCVL  = %1d\nkLong = %1d\n',kYOU,kCVL,kLONG)

% Example 2
alfa = 0.0;
beta = 0.9;
kYOU = nEigs(U,alfa,beta);
kCVL = nEigsCVL(U,alfa,beta);
ev = real(eig(U));
kLONG = sum(alfa-tol<=ev & ev<=beta+tol);
fprintf('\nm = %2d   n = %2d  alfa = %6.3f  beta = %6.3f\n',m,n,alfa,beta)
fprintf('kYOU  = %1d\nkCVL  = %1d\nkLong = %1d\n',kYOU,kCVL,kLONG)

m = 100;
n = 20;
[U,R] = qr(randn(m,n));
tol = 10^(-12);

% Example 3
alfa = -0.999;
beta =  0.999;
kYOU = nEigs(U,alfa,beta);
kCVL = nEigsCVL(U,alfa,beta);
ev = real(eig(U));
kLONG = sum(alfa-tol<=ev & ev<=beta+tol);
fprintf('\nm = %2d   n = %2d  alfa = %6.3f  beta = %6.3f\n',m,n,alfa,beta)
fprintf('kYOU  = %1d\nkCVL  = %1d\nkLong = %1d\n',kYOU,kCVL,kLONG)




% Example 4
alfa = -1.0;
beta = -1.0;
kYOU = nEigs(U,alfa,beta);
kCVL = nEigsCVL(U,alfa,beta);
ev = real(eig(U));
kLONG = sum(alfa-tol<=ev & ev<=beta+tol);
fprintf('\nm = %2d   n = %2d  alfa = %6.3f  beta = %6.3f\n',m,n,alfa,beta)
fprintf('kYOU  = %1d\nkCVL  = %1d\nkLong = %1d\n',kYOU,kCVL,kLONG)


% Example 5
alfa = 1.0;
beta = 1.0;
kYOU = nEigs(U,alfa,beta);
kCVL = nEigsCVL(U,alfa,beta);
ev = real(eig(U));
kLONG = sum(alfa-tol<=ev & ev<=beta+tol);
fprintf('\nm = %2d   n = %2d  alfa = %6.3f  beta = %6.3f\n',m,n,alfa,beta)
fprintf('kYOU  = %1d\nkCVL  = %1d\nkLong = %1d\n',kYOU,kCVL,kLONG)




  function k = nEigsCVL(U,alfa,beta)
% U is an nxn real orthogonal matrix.
% If -1 < alfa < beta < 1, then k estimates the number of U's eigenvalues that 
%          have real parts in the interval [alfa,beta].
% If alfa = beta = 1, then k estimates the number of U's eigenvalues that
%     are equal to one.
% If alfa = beta = -1, then k estimates the number of U's eigenvalues that
%     are equal to minus one.
% All other alfa-beta combinations are illegal, e.g., alfa = 0, beta = 1

[n,n] = size(U);
tol = 10^-12;
if alfa==1 && beta==1
    % In this case, k = n - rank(U-I). Use QR with col pivoting.
    [Q,R,P] = qr(U-eye(n,n));
    k = n-sum(abs(diag(R))>=tol);
    return
end
if alfa==-1 && beta==-1
    % In this case, k = n - rank(U+I). Use QR with col pivoting.
    [Q,R,P] = qr(U+eye(n,n));
    k = n-sum(abs(diag(R))>=tol);
    return
end
T = hess((U+U')/2);
% Where are the zero subdiagonal?
kdx = find(abs(diag(T,-1))<=tol);
% H(kdx(i)+1,kdx(i)) is small for i=1:length(kdx)
First = 1;
k = 0;
for j=1:length(kdx)+1
    if j>length(kdx)
        Last = n;
    else
        Last = kdx(j);
    end
    ran = First:Last;
    k = k+nLess(T(ran,ran),beta)-nLess(T(ran,ran),alfa);
    First=Last+1;
end