ulv_qrit.m

UTV工具包提供46个Matlab函数

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.')
end
if (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
  end
end

%-----------------------------------------------------------------------
% End of function ulv_qrit
%-----------------------------------------------------------------------