function [Q,D] = SpecialSchur(X,B)
% X is an n-by-p matrix, 1<=p<n
% B is a symmetric p-by-p matrix
% A = I + X*B*X'
% Q'*A*Q = D is the schur decomposition of A.

[n,p] = size(X);
[Q,R] = qr(X);
R1 = R(1:p,1:p);
Btilde = R1*B*R1';        % R*B*R' is O(pn^2)    R1*B*R1' is O(p^3)
% Note Atilde = Q'*(I + XBX')Q 
%             = I + (Q'X)B(Q'X)'  =  eye(p,p)+Btilde  0; 0 eye(n-p,n-p)]

[U1,D1] = schur(eye(p,p)+R1*B*R1');
% Q = Q*[U1 0 ; 0 eye(n-p)]
Q(:,1:p) = Q(:,1:p)*U1;
D = [D1 zeros(p,n-p) ; zeros(n-p,p) eye(n-p,n-p)];

