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