function [L,V,U] = ulv_qrit(p,num_ref,L,V,U)%  ulv_qrit --> Refinement of L in the ULV decomposition via QR-iterations.%%  <Synopsis>%    [L,V,U] = ulv_qrit(p,num_ref,L,V,U)%    [L,V]   = ulv_qrit(p,num_ref,L,V)%    [L]     = ulv_qrit(p,num_ref,L)%%  <Description>%    Given the ULV decomposition U*L*V' with numerical rank p, the%    function refines the rank-revealing decomposition via num_ref%    steps of block QR iterations.%%  <Algorithm>%    Refinement is identical to block QR iteration, in which the%    off-diagonal block of the lower triangular matrix L is "flipped"%    to the (1,2)-position and then back again.%%  <See Also>%    ulv_ref --> Refine one row of L in the ULV decomposition.%  <References>%  [1] R. Mathias and G.W. Stewart, "A Block QR Algorithm and the Singular%      Value Decomposition", Lin. Alg. Appl., 182 (1993), pp. 91--100.%%  <Revision>%    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 number of input and output arguments.vflag = 1;uflag = 1;if (nargout == 1 & nargin == 3)  vflag = 0;  uflag = 0;elseif (nargout == 2 & nargin == 4)  uflag = 0;elseif ~(nargout == 3 & nargin == 5)  error('Wrong number of input or output arguments.')end% Initialize.[n,n] = size(L);if (p ~= abs(round(p))) | (p > n)  error('Requires p to be an integer between 0 and n.')endif (num_ref ~= abs(round(num_ref)))  error('Requires positive integer value for num_ref.')end% Refinement loop.for (i = 1:num_ref)  % Flip the off-diagonal block of L using Givens rotations on the left.  for (r = p+1:n)    for (i = p:-1:1)      % Apply rotation to L on the left.      [c,s,L(i,i)] = gen_giv(L(i,i),L(r,i));      L(r,i) = 0;                               % Eliminate L(r,i).      [L(i,1:i-1),L(r,1:i-1)] = app_giv(L(i,1:i-1),L(r,1:i-1),c,s);      [L(i,p+1:r),L(r,p+1:r)] = app_giv(L(i,p+1:r),L(r,p+1:r),c,s);      % Apply rotation to U on the right.      if (uflag)        [U(:,i),U(:,r)] = app_giv(U(:,i),U(:,r),c,s);      end    end  end  % Restore lower triangular form with Givens rotations on the right.  for (r = n:-1:p+1)    for (i = 1:p)      % Apply rotation to L on the right.      [c,s,L(i,i)] = gen_giv(L(i,i),L(i,r));      L(i,r) = 0;                               % Eliminate L(i,r).      [L(i+1:p,i),L(i+1:p,r)] = app_giv(L(i+1:p,i),L(i+1:p,r),c,s);      [L(r:n,i),L(r:n,r)]     = app_giv(L(r:n,i),L(r:n,r),c,s);      % Apply rotation to V on the right.      if (vflag)        [V(1:n,i),V(1:n,r)] = app_giv(V(1:n,i),V(1:n,r),c,s);      end    end  endend%-----------------------------------------------------------------------% End of function ulv_qrit%-----------------------------------------------------------------------