lanc_b.m

这是在网上下的一个东东

function [U,B_k,V] = lanc_b(A,p,k,reorth) 
%LANC_B Lanczos bidiagonalization. 
% 
% B_k = lanc_b(A,p,k,reorth) 
% [U,B_k,V] = lanc_b(A,p,k,reorth) 
% 
% Performs k steps of the Lanczos bidiagonalization process with 
% starting vector p, producing a lower bidiagonal matrix 
%           [b_11               ] 
%           [b_21 b_22          ] 
%     B_k = [     b_32 .        ] 
%           [          . b_kk   ] 
%           [            b_k+1,k] 
% such that 
%     A*V = U*B_k , 
% where U and V consist of the left and right Lanczos vectors. 
% 
% Reorthogonalization is controlled by means of reorth: 
%    reorth = 0 : no reorthogonalization, 
%    reorth = 1 : reorthogonalization by means of MGS, 
%    reorth = 2 : Householder-reorthogonalization. 
% No reorthogonalization is assumed if reorth is not specified. 
 
% Reference: G. H. Golub & C. F. Van Loan, "Matrix Computations", 
% 3. Ed., Johns Hopkins, 1996.  Section 9.3.4. 
% Referred to as "bidiag1" by Paige and Saunders. 
 
% Per Christian Hansen, IMM, April 8, 2001. 
 
% Initialization. 
if (k<1), error('Number of steps k must be positive'), end 
if (nargin < 4), reorth = 0; end 
if (reorth < 0 | reorth > 2), error('Illegal reorth'), end 
if (nargout==2), error('Not enough output arguments'), end 
[m,n] = size(A); 
B_k = sparse(k+1,k); 
if (nargout>1 | reorth==1) 
  U = zeros(m,k); V = zeros(n,k); UV = 1; 
else 
  UV = 0; 
end 
if (reorth==2) 
  if (k>=n), error('No. of iterations must satisfy k < n'), end 
  HHU = zeros(m,k); HHV = zeros(n,k); 
  HHalpha = zeros(1,k); HHbeta = HHalpha; 
end 
 
% Prepare for Lanczos iteration. 
v = zeros(n,1); 
beta = norm(p); 
if (beta==0), error('Starting vector must be nonzero'), end 
if (reorth==2) 
  [beta,HHbeta(1),HHU(:,1)] = gen_hh(p); 
end 
u = p/beta; 
if (UV), U(:,1) = u; end 
 
% Perform Lanczos bidiagonalization with/without reorthogonalization. 
for i=1:k 
 
  r = (u'*A)' - beta*v;   % A'*u 
  if (reorth==0) 
    alpha = norm(r); v = r/alpha; 
  elseif (reorth==1) 
    for j=1:i-1, r = r - (V(:,j)'*r)*V(:,j); end 
    alpha = norm(r); v = r/alpha; 
  else 
    for j=1:i-1 
      r(j:n) = app_hh(r(j:n),HHalpha(j),HHV(j:n,j)); 
    end 
    [alpha,HHalpha(i),HHV(i:n,i)] = gen_hh(r(i:n)); 
    v = zeros(n,1); v(i) = 1; 
    for j=i:-1:1 
      v(j:n) = app_hh(v(j:n),HHalpha(j),HHV(j:n,j)); 
    end 
  end 
  B_k(i,i) = alpha; if (UV), V(:,i) = v; end 
 
  p = A*v - alpha*u; 
  if (reorth==0) 
    beta = norm(p); u = p/beta; 
  elseif (reorth==1) 
    for j=1:i, p = p - (U(:,j)'*p)*U(:,j); end 
    beta = norm(p); u = p/beta; 
  else 
    for j=1:i 
      p(j:m) = app_hh(p(j:m),HHbeta(j),HHU(j:m,j)); 
    end 
    [beta,HHbeta(i+1),HHU(i+1:m,i+1)] = gen_hh(p(i+1:m)); 
    u = zeros(m,1); u(i+1) = 1; 
    for j=i+1:-1:1 
      u(j:m) = app_hh(u(j:m),HHbeta(j),HHU(j:m,j)); 
    end 
  end 
  B_k(i+1,i) = beta; if (UV), U(:,i+1) = u; end 
 
end 
 
if (nargout==1), U = B_k; end