bidiag.m

这是在网上下的一个东东

function [U,B,V] = bidiag(A) 
%BIDIAG Bidiagonalization of an m-times-n matrix with m >= n. 
% 
% B = bidiag(A) 
% [U,B,V] = bidiag(A) 
%  
% Computes the bidiagonalization of the m-times-n matrix A with m >= n: 
%     A = U*B*V' , 
% where B is an upper bidiagonal n-times-n matrix, and U and V have 
% orthogonal columns.  The matrix B is stored as a sparse matrix. 
 
% Reference: L. Elden, "Algorithms for regularization of ill- 
% conditioned least-squares problems", BIT 17 (1977), 134-145. 
 
% Per Christian Hansen, IMM, Sept. 13, 2001.
 
% Initialization. 
[m,n] = size(A); 
if (m < n), error('Illegal dimensions of A'), end 
B = sparse(n,n); 
if (nargout> 1), U = [eye(n);zeros(m-n,n)]; betaU = zeros(n,1); end 
if (nargout==3), V = eye(n); betaV = zeros(n,1); end 
 
% Bidiagonalization; save Householder quantities. 
if (m > n), k_last = n; else k_last = n-1; end 
for k=1:k_last 
 
  [B(k,k),beta,A(k:m,k)] = gen_hh(A(k:m,k)); 
  if (k < n), A(k:m,k+1:n) = app_hh(A(k:m,k+1:n),beta,A(k:m,k)); end 
  if (nargout>1), betaU(k) = beta; end 
 
  if (k < n-1) 
    [B(k,k+1),beta,v] = gen_hh(A(k,k+1:n).'); A(k,k+1:n) = v.'; 
    A(k+1:m,k+1:n) = app_hh(A(k+1:m,k+1:n)',beta,A(k,k+1:n)')'; 
    if (nargout==3), betaV(k) = beta;, end 
  elseif (k == n-1) 
    B(n-1,n) = A(n-1,n); 
  end 
 
end 
 
% Save bottom element if A is square. 
if (k_last < n), B(n,n) = A(n,n); end 
 
% Compute U if wanted. 
if (nargout>1) 
  for k=k_last:-1:1 
    U(k:m,k:n) = app_hh(U(k:m,k:n),betaU(k),A(k:m,k)); 
  end 
end 
 
% Compute V if wanted. 
if (nargout==3) 
  for k=n-2:-1:1 
    V(k+1:n,k:n) = app_hh(V(k+1:n,k:n),betaV(k),A(k,k+1:n)'); 
  end 
end 
 
if (nargout < 2), U = B; end