function [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank,fixed_rank)%  hrrqr --> Chan/Foster high-rank-revealing RRQR algorithm.%%  <Synopsis>%    [p,R,Pi,Q,W,vec] = hrrqr(A)%    [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank)%    [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank,fixed_rank)%%  <Description>%    Computes a rank-revealing RRQR decomposition of an m-by-n matrix A%    (m >= n) with numerical rank p close to n. The n-by-n matrix R is%    upper triangular and will reveal the numerical rank p of A. Thus,%    the norm of the (2,2) block of R is of the order sigma_(p+1).%%  <Input Parameters>%    1. A          --> m-by-n matrix (m >= n);%    2. tol_rank   --> rank decision tolerance;%    3. fixed_rank --> deflate to the fixed rank given by fixed_rank instead%                      of using the rank decision tolerance;%%    Defaults: tol_rank = sqrt(n)*norm(A,1)*eps;%%  <Output parameters>%    1.   p        --> numerical rank of A;%    2-4. R, Pi, Q --> the RRQR factors so that A*Pi = Q*R;%    5.   W        --> an n-by-p matrix whose columns span an%                      approximation to the null space of A;%    6.   vec      --> a 5-by-1 vector with:%         vec(1) = upper bound of norm(R(1:p,p+1:n)),%         vec(2) = estimate of pth singular value,%         vec(3) = estimate of (p+1)th singular value,%         vec(4) = a posteriori upper bound of num. nullspace angle,%         vec(5) = a posteriori upper bound of num. range angle.%%  <Algorithm>%    The rectangular matrix A is preprocessed by a QR factorization, A = Q*R.%    Then deflation steps based on the generalized LINPACK condition%    estimator are employed to produce a rank-revealing decomposition.%%  <See Also>%    lrrqr --> Chan-Hansen low-rank-revealing RRQR algorithm.%  <References>%  [1] T.F. Chan, "Rank Revealing QR Factorizations", Lin.%      Alg. Appl., 88/89 (1987), pp. 67--82.%%  [2] L. Foster, "Rank and Null Space Calculations Using Matrix Decomposition%      Without Column Interchanges", Lin. Alg. Appl., 74 (1986), pp. 47--71.%%  <Revision>%    Christian H. Bischof, Argonne National Laboratory%    Ricardo D. Fierro, California State University San Marcos%    Per Christian Hansen, IMM, Technical University of Denmark%    Peter S.K. Hansen, IMM, Technical University of Denmark%%    Last revised: June 22, 1999%-----------------------------------------------------------------------% Check the required input arguments.if (nargin < 1)  error('Not enough input arguments.')end[m,n] = size(A);if (m*n == 0)  error('Empty input matrix A not allowed.')elseif (m < n)  error('The system is underdetermined.')end% Check the optional input arguments, and set defaults.if (nargin == 1)  tol_rank   = sqrt(n)*norm(A,1)*eps;  fixed_rank = 0;elseif (nargin == 2)  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end  fixed_rank = 0;elseif (nargin == 3)  if isempty(fixed_rank)    if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end    fixed_rank = 0;  else    tol_rank = realmax;    if (fixed_rank ~= abs(round(fixed_rank))) | (fixed_rank > n)      error('Requires fixed_rank to be an integer between 0 and n.')    end  endendif (tol_rank ~= abs(tol_rank))  error('Requires positive value for tol_rank.')end% Check the number of output arguments.piflag  = 1;qflag   = 1;wflag   = 1;vecflag = 1;if (nargout <= 2)  piflag  = 0; Pi = [];  qflag   = 0; Q  = [];  wflag   = 0; W  = [];  vecflag = 0;elseif (nargout == 3)  qflag   = 0; Q = [];  wflag   = 0; W = [];  vecflag = 0;elseif (nargout == 4)  wflag   = 0; W = [];  vecflag = 0;elseif (nargout == 5)  vecflag = 0;end% Compute initial skinny QR factorization A = Q*R.if (qflag)  [Q,R] = qr(A,0);else  R = triu(qr(A));  R = R(1:n,1:n);endif (piflag)  Pi = eye(n);endif (wflag)  W = zeros(n,0);end% Rank-revealing procedure.% Estimate of the n'th singular value and the corresponding right% singular vector via the generalized LINPACK condition estimator.[smin,v] = ccvl(R);smin_p_plus_1 = 0;                             % No (n+1)th singular value.p = n;                                         % Init. loop to full rank n.while ((smin < tol_rank) & (p > fixed_rank))  % Update the matrix W.  if (wflag)    W = [[v;zeros(n-p,1)],W];  end  % Find the element in v with greatest index in absolute value.  [vmax,k] = max(abs(v)); k = max(k);  % If necessary, generate the permutation that brings the pivot  % element of v to the last position, apply the permutation to W,  % Pi, and R, compute a new QR factorization of R(1:i,1:i), and  % update Q and R.  if (k < p)    perm = [k+1:p,k];    R(:,k:p) = R(:,perm);    if (piflag)      Pi(:,k:p) = Pi(:,perm);    end    if (wflag)      W(k:p,:) = W(perm,:);    end    for (j = k:p-1)      [c,s,R(j,j)] = gen_giv(R(j,j),R(j+1,j)); R(j+1,j) = 0;      [R(j,j+1:n),R(j+1,j+1:n)] = app_giv(R(j,j+1:n),R(j+1,j+1:n),c,s);      if (qflag)        [Q(:,j),Q(:,j+1)] = app_giv(Q(:,j),Q(:,j+1),c,s);      end    end  end  % New rank estimate after the problem has been deflated.  p = p - 1;  smin_p_plus_1 = smin;  % Estimate of the p'th singular value and the corresponding right  % singular vector via the generalized LINPACK condition estimator.  if (p > 0)    [smin,v] = ccvl(R(1:p,1:p));  else    smin = 0;                                  % No 0th singular value.  endendif (wflag)  W = Pi*W;end% Estimates that describe the quality of the decomposition.if (vecflag)  vec = zeros(5,1);  if (p > 0)    vec(1) = sqrt(n-p)*norm(R(1:p,p+1:n),1);    vec(2) = smin;    vec(3) = smin_p_plus_1;    vec(4) = (smin_p_plus_1^2)/(smin^2 - smin_p_plus_1^2);    vec(5) = smin_p_plus_1/smin;  else    vec(3) = smin_p_plus_1;  endend%-----------------------------------------------------------------------% End of function hrrqr%-----------------------------------------------------------------------