function [U,T,r,berr] = rrtrid(A,tol)
%RRTRID   Rank-revealing tridiagonalization.
%   T = RRTRID(A,tol) is the tridiagonal form of the 
%   (skew-)symmetric/Hermitian matrix A with a trailing zero diagonal 
%   block in case of numerical rank deficiency.  
%   By default, tol = max(size(A)) * eps(norm(A)) as in 'rank'.
%   [U,T] = RRTRID(A) produces a unitary matrix U and a tridiagonal
%   matrix T so that A = U*T*U' and T'*T = EYE(SIZE(T)).
%   [U,T,r] = RRTRID(A) additionally returns the numerical rank r
%   of A with respect to tol as computed by RRQR.
%   [U,T,r,berr] = RRTRID(A) also returns the relative backward error 
%   introduced by setting to zero all matrix entries that are neglected 
%   due to RRQR's rank determination.  The erros are measured in the
%   2-norm.  Use this for diagnostics only. 
%   Note that on output, all these elements are explicitly set to zero.  

% Peter Benner, 2007-11-15

[n,m] = size(A);
if n~=m,  error('A must be square.'),  end
issymm = 0;
isskew = 0;
if (norm(A-A',1) == 0),
    issymm = 1;
elseif (norm(A+A',1) == 0),
    isskew = 1;
else
    error('A must be symmetric (Hermitian) or skew-symmetric(-Hermitian).')
end
if ~isreal(A),
    error('Currently, only real arithmetic suppoerted.')
end
nrmA = norm(A);
if (nargin < 2) || (tol < 0),
    tol = n * eps(nrmA);
end
[U,R,p,r] = rrqr(A,tol);
T = R(1:r,:)*U(p,:);
if nargout > 3,
    berr = (norm(R(r+1:end,:)*U(p,:)) + norm(T(:,r+1:end)))/nrmA;
end
T = T(1:r,1:r);
if issymm,
    T = (T + T')/2;
else
    T = (T - T')/2;
end  
[P,T] = hess(T);
%hess does not respect skew structures.
if isskew,
    T = T - triu(T,2) - diag(diag(T));
end
T = [T zeros(r,n-r); zeros(n-r,n)];
if nargout > 1,
    U(:,1:r) = U(:,1:r)*P;
else
    U = T;
end