function P4Grade()
clc
% Small example for correctness
n = 20; r = 5; X = randn(n,r); Y = randn(n,r); order = 'descend';
[U,T]  = specialSchur(X,Y,order);
check(U,T,X,Y,order)
fprintf('\n\n')
% Small example for correctness
n = 20; r = 5; X = randn(n,r); Y = randn(n,r); order = 'ascend';
[U,T]  = specialSchur(X,Y,order);
check(U,T,X,Y,order)
fprintf('\n\n')
% Large example for efficiency...
n = 1000; r = 5; X = randn(n,r); Y = randn(n,r); order = 'descend';
tic,[U,T1] = specialSchurCVL1(X,Y,order); tCVL1 = toc; % check(U,T3,X,Y,order)
tic,[U,T2] = specialSchurCVL2(X,Y,order); tCVL2 = toc; % check(U,T3,X,Y,order)
tic,[U,T3] = specialSchurCVL3(X,Y,order); tCVL3 = toc; % check(U,T3,X,Y,order)
tic,[U,T]  = specialSchur(X,Y,order); tYOU = toc; check(U,T,X,Y,order)
disp(' ')
fprintf('tYOU/tCVL1 = %6.3f\n',tYOU/tCVL1)
fprintf('tYOU/tCVL2 = %6.3f\n',tYOU/tCVL2)
fprintf('tYOU/tCVL3 = %6.3f\n',tYOU/tCVL3)



function [U,T] = specialSchurCVL1(X,Y,order)
% X and Y are n-by-r matrices with r<n/2. 
% U is orthogonal and U'*(X*Y')*U = T is the real Schur decomposition of XY'.
% If order = 'ascend' then the eigenvalues along the block diagonal of T are arranged
%       in order of increasing real part.
% If order = 'descend' then the eigenvalues along the block diagonal of T are arranged
%       in order of decreasing real part.
%

% This implementation disregards the low rank structure...
[U,T] = schur(X*Y');
[U,T] = reorder(U,T,order);


function [U,T] = specialSchurCVL2(X,Y,order)
% X and Y are n-by-r matrices with r<n/2. 
% U is orthogonal and U'*(X*Y')*U = T is the real Schur decomposition of XY'.
% If order = 'ascend' then the eigenvalues along the block diagonal of T are arranged
%       in order of increasing real part.
% If order = 'descend' then the eigenvalues along the block diagonal of T are arranged
%       in order of decreasing real part.

[n,r] = size(X);
[U,R] = qr([X,Y]);
%                    R11  R12
% Since [U'X U'Y] =   0   R22
%                     0    0
% where
R11 = R(1:r,1:r); R12 = R(1:r,r+1:2*r); R22 = R(r+1:2*r,r+1:2*r);
% It follows that
%                        R11*R12'   R11*R22'  0
%           U'(XY')U  =     0          0      0
%                           0          0      0
% Compute the rxr real Schur decomposition...
[Q11,T11] = schur(R11*R12');
% If
U(:,1:r) = U(:,1:r)*Q11;
T12 = Q11'*R11*R22';
% then
%                        T11  T12  0
%         U'*(XY')U =    0     0   0  = T
%                        0     0   0
T = zeros(n,n); T(1:r,1:r) = T11; T(1:r,r+1:2*r) = T12;
[U,T] = reorder(U,T,order);


function [U,T] = specialSchurCVL3(X,Y,order)
% X and Y are n-by-r matrices with r<n/2. 
% U is orthogonal and U'*(X*Y')*U = T is the real Schur decomposition of XY'.
% If order = 'ascend' then the eigenvalues along the block diagonal of T are arranged
%       in order of increasing real part.
% If order = 'descend' then the eigenvalues along the block diagonal of T are arranged
%       in order of decreasing real part.
%
%  This is just like specialSchurCVL2, except that instead of calling
%  ordschur on an nxn problem, we realize the necessary reordering via a
%  block-level permutation and a pair of O(r) calls to ordschur.

[n,r] = size(X);
[U,R] = qr([X,Y]);
%                    R11  R12
% Since [U'X U'Y] =   0   R22
%                     0    0
% where
R11 = R(1:r,1:r); R12 = R(1:r,r+1:2*r); R22 = R(r+1:2*r,r+1:2*r);
% It follows that
%                        R11*R12'   R11*R22'  0
%           U'(XY')U  =     0          0      0
%                           0          0      0
% Compute the rxr real Schur decomposition...
[Q11,T11] = schur(R11*R12');
% Reorder...
[Q11,T11] = reorder(Q11,T11,order);
% If
U(:,1:r) = U(:,1:r)*Q11;
T12 = Q11'*R11*R22';
% then
%                        T11  T12  0
%         U'*(XY')U =    0     0   0  = T
%                        0     0   0
T = zeros(n,n); T(1:r,1:r) = T11; T(1:r,r+1:2*r) = T12;



% Some of T11's eigs have to be sent the lower right corner of T.
% "Re-block"  T:
%        T11 T12  T13  0    T11:  r1-x-r1  T12:  r1-x-(r-r1)   T13: r1-x-r
%    T =  0  T22  T23  0    T22: (r-r1)-x-(r-r1)     T23: (r-r1)-x-r
%         0   0    0   0
%         0   0    0   0
% And move the T22 eigs to the bottom. Towards this end,
% interchange block rows and columns 2 and 4 and this gives

e = ordeig(T11);
tol = .000001;
if strcmp(order,'ascend')
    r1 = sum(real(e)<=0 & abs(e)>=tol);
else
    r1 = sum(real(e)>=0 & abs(e)>=tol);
end
idx = [1:r1 (2*r+1):n (r+1):(2*r)  (r1+1):r];
T = T(idx,idx);
U = U(:,idx);
%                     T11   0   T13  T12
%           T <---     0    0    0    0
%                      0    0    0    0
%                      0    0   T23  T22
% and 
%        U'(XY')U = T
%
% We apply Schur to
%                           0    0
%                          T23  T22
m = r + r-r1;
ran = (n-m+1):n;
[Q22,T22] = schur(T(ran,ran));
[Q22,T22] = reorder(Q22,T22,order);
% and update...
U(:,ran) = U(:,ran)*Q22;
T(ran,ran) = T22;
T(1:r1,ran) = T(1:r1,ran)*Q22;


   
function [U1,T1] = reorder(U,T,order)
% Does an ordschur update of U and T according to whether
% order = 'ascend' or 'descend'
% We need a definition of zero in order to decide what kind of
% eigenvalue we have.
tol = .000001;
% The eigenvalues of T...
d = ordeig(T);
% Find out where the zero eigs are and set them to zero..
zero_idx    = find(abs(d)<=tol);
%d(zero_idx) = 0;
% Find out where the left half plane and right half plane eigs are...
left_idx    = find(real(d)<=0 & abs(d)>tol);
right_idx   = find(real(d)> 0 & abs(d)>tol);
% Set up the select vector for ordschur...
n = length(d);
select = zeros(n,1);
if strcmp(order,'ascend')
   select(left_idx)  = 3;
   select(zero_idx)  = 2;
   select(right_idx) = 1;
else
   select(left_idx)  = 1;
   select(zero_idx)  = 2;
   select(right_idx) = 3;
end
% Reorder...
[U1,T1] = ordschur(U,T,select);


  function check(U,T,X,Y,order)
% Reports on the output of specialSchur
[n,r] = size(X);
tol = eps*norm(X,'fro')*norm(Y,'fro');
tol = .000001;
% U is orthogonal...
errU = norm(U'*U - eye(n,n));
% T is similar to XY' via U
errT = norm(U'*(X*Y')*U - T);
% The eigenvalues as they appear along T's diagonal...
e = ordeig(T);
% Where are the right half plane, zero, and left half plane eigenvalues?
nPos    = sum( real(e)>=0 & abs(e)>=tol);
idx_Pos = find(real(e)>=0 & abs(e)>=tol); 
nZer    = sum(abs(e)<tol);
idx_Zer = find(abs(e)<tol);
nNeg    = sum( real(e)<0 & abs(e)>=tol);
idx_Neg = find(real(e)<0 & abs(e)>=tol); 
if strcmp(order,'descend')
    % right-to left...
    PosOK = (length(idx_Pos)==nPos) && norm(sort(idx_Pos) - (1:nPos)')==0;
    ZerOK = (length(idx_Zer)==nZer) && norm(sort(idx_Zer) - (nPos+1:nPos+nZer)')==0;
    NegOK = (length(idx_Neg)==nNeg) && norm(sort(idx_Neg) - (n-nNeg+1:n)')==0;
else
    % left-to-right...
    NegOK = (length(idx_Neg)==nNeg) && norm(sort(idx_Neg) - (1:nNeg)')==0;
    ZerOK = (length(idx_Zer)==nZer) && norm(sort(idx_Zer) - (nNeg+1:nNeg+nZer)')==0;
    PosOK = (length(idx_Pos)==nPos) && norm(sort(idx_Pos) - (n-nPos+1:n)')==0;
end
fprintf('n = %4d\nr = %4d\norder = %s\n||U''U - I||       = %10.3e\n||U''(XY'')U - T || = %10.3e\n',...
              n,r,order,errU,errT)
if PosOK
    fprintf('nPos = %3d  PosOK = TRUE\n',nPos)
else
    fprintf('nPos = %3d  PosOK = FALSE\n',nPos)
end
if ZerOK
    fprintf('nZer = %3d  ZerOK = TRUE\n',nZer)
else
    fprintf('nZer = %3d  ZerOK = FALSE\n',nZer)
end
if NegOK
    fprintf('nNeg = %3d  NegOK = TRUE\n',nNeg)
else
    fprintf('nNeg = %3d  NegOK = FALSE\n',nNeg)
end


  