ccvl.m

UTV工具包提供46个Matlab函数

function [smin,vmin] = ccvl(R)
%  ccvl --> Singular value/vector estimates via condition estimation.
%
%  <Synopsis>
%    [smin,vmin] = ccvl(R)
%
%  <Description>
%    Compute estimates smin and vmin of the smallest singular value
%    and corresponding right singular vector of the upper triangular
%    matrix R.
%
%  <Algorithm>
%    The function is based on the generalized LINPACK condition
%    number estimator.
%
%  <See Also>
%    inviter --> Singular value/vector estimates via inverse iteration.

%  <References>
%  [1] A.K. Cline, A.R. Conn & C.F. Van Loan, "Generalizing the LINPACK
%      Condition Estimator"; in J.P. Hennart (Ed.), "Numerical Analysis",
%      Lecture Notes in Mathematics, Vol. 909, pp. 73-83, Springer, (1982).
%
%  <Revision>
%    Christian H. Bishof, Argonne National Laboratory
%    C.-T. Pan, NIU
%    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
%-----------------------------------------------------------------------

% Initialize.
[n,n] = size(R);
if (n==0)
  smin = [];
  vmin = [];
  return
end

% Matrix is singular.
% if (any(diag(R)==0))
%   smin = 0;
%   vmin = zeros(n,1);
%   vmin(min(find(diag(R)==0))) = 1;
%   return
% end

for i=1:n, if R(i,i)==0, R(i,i) = eps; end, end

wsave = warning;
warning('off');

v = zeros(n,1);

% First element is treated special.
v(1)  = 1/R(1,1);
vnorm = abs(v(1));
p     = v(1)*R(1,:);

% Process rows of matrix one by one.
for (i = 2:n-1)
  u     = R(i,i+1:n);
  utp   = u*p(i+1:n)';
  gamma = R(i,i);
  xi    = p(i);
  phi   = 1 + norm(u)^2;
  pnorm = norm(p(i+1:n));
  alpha = xi*phi - gamma*utp;

  if (alpha == 0)
    beta = gamma^2*(vnorm^2 + pnorm^2) - (xi^2+1)*phi;
    if (beta > 0)
      s = 1;
      c = 0;
    else
      s = 0;
      c = 1;
    end
  else
    beta = gamma^2*(vnorm^2 + pnorm^2) + (xi^2-1)*phi - 2*xi*gamma*utp;
    eta  = beta/pythag(beta,2*alpha);       % eta is cos(2a).
    s    = -sign(alpha)*sqrt((1+eta/2));
    c    = sqrt((1-eta)/2);
  end

  v(1:i)   = [s*v(1:i-1);(c-s*xi)/gamma];
  vnorm    = pythag(s*vnorm,v(i));
  p(i+1:n) = s*p(i+1:n) + v(i)*u;
end

% Last step is the same as in incremental condition estimation.
alpha = p(n);
gamma = R(n,n);
if (alpha == 0)
  beta = gamma^2*vnorm^2 - 1;
  if (beta > 0)
    s = 1;
    c = 0;
    lambda_max = beta + 1;
  else
    s = 0;
    c = 1;
    lambda_max = 1;
  end
else
  beta = gamma^2*vnorm^2 + alpha^2 - 1;
  eta  = beta/pythag(beta,2*alpha);
  lambda_max = 0.5*(beta + pythag(beta,2*alpha)) + 1;
  s = -sign(alpha)*sqrt((1+eta)/2);
  c = sqrt((1-eta)/2);
end
v     = [s*v(1:n-1);(c-s*alpha)/gamma];
vnorm = sqrt(lambda_max)/abs(gamma);

% Compute and normalize right nullvector.
vmin = R\(v/vnorm);
smin = 1/norm(vmin);
vmin = smin*vmin;

warning(wsave);

%-----------------------------------------------------------------------
% End of function ccvl
%-----------------------------------------------------------------------

function x = pythag(y,z)
%  pythag --> Pythagoras equation.
%
%  <Synopsis>
%    x = pythag(y,z)
%
%  <Description>
%    Returns sqrt(y*y + z*z) but is careful to scale to avoid overflow.

%  <Revision>
%    Christian H. Bishof, Argonne National Laboratory
%
%    Revised: March 31, 1989
%-----------------------------------------------------------------------

rmax = max(abs([y;z]));
if (rmax == 0)
  x = 0;
else
  x = rmax*sqrt((y/rmax)^2 + (z/rmax)^2);
end

%-----------------------------------------------------------------------
% End of function pythag
%-----------------------------------------------------------------------