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