www.pudn.com > RobustSF.zip > lanpro.m, change:2000-05-10,size:14762b

function [Q_k,T_k,r,anorm,ierr,work] = lanpro(A,nin,kmax,r,options,...
Q_k,T_k,anorm)

%LANPRO   Lanczos tridiagonalization with partial reorthogonalization
%   LANPRO computes the Lanczos tridiagonalization of a real symmetric
%   matrix using the symmetric Lanczos algorithm with partial
%   reorthogonalization.
%
%   [Q_K,T_K,R,ANORM,IERR,WORK] = LANPRO(A,K,R0,OPTIONS,Q_old,T_old)
%   [Q_K,T_K,R,ANORM,IERR,WORK] = LANPRO('Afun',N,K,R0,OPTIONS,Q_old,T_old)
%
%   Computes K steps of the Lanczos algorithm with starting vector R0,
%   and returns the K x K tridiagonal T_K, the N x K matrix Q_K
%   with semiorthonormal columns and the residual vector R such that
%
%        A*Q_K = Q_K*T_K + R .
%
%   Partial reorthogonalization is used to keep the columns of Q_K
%   semiorthogonal:
%        MAX(DIAG((eye(k) - Q_K'*Q_K))) <= OPTIONS.delta.
%
%
%   The first input argument is either a real symmetric matrix, a struct with
%   components A.L and A.U or a string containing the name of an M-file which
%   applies a linear operator to the columns of a given matrix.  In the latter
%   case, the second input argument must be N, the order of the problem.
%
%   If A is a struct with components A.L and A.U, such that
%   L*U = (A - sigma*I), a shift-and-invert Lanczos iteration is performed
%
%   The OPTIONS structure is used to control the reorthogonalization:
%     OPTIONS.delta:  Desired level of orthogonality
%                     (default = sqrt(eps/K)).
%     OPTIONS.eta  :  Level of orthogonality after reorthogonalization
%                     (default = eps^(3/4)/sqrt(K)).
%     OPTIONS.cgs  :  Flag for switching between different reorthogonalization
%                     algorithms:
%                      0 = iterated modified Gram-Schmidt  (default)
%                      1 = iterated classical Gram-Schmidt
%     OPTIONS.elr  :  If OPTIONS.elr = 1 (default) then extended local
%                     reorthogonalization is enforced.
%     OPTIONS.Y    :  The lanczos vectors are reorthogonalized against
%                     the columns of the matrix OPTIONS.Y.
%
%   If both R0, Q_old and T_old are provided, they must contain
%   a partial Lanczos tridiagonalization of A on the form
%
%        A Q_old = Q_old T_old + R0 .
%
%   In this case the factorization is extended to dimension K x K by
%   continuing the Lanczos algorithm with R0 as starting vector.
%
%   On exit ANORM contains an approximation to ||A||_2.
%     IERR = 0  :  K steps were performed succesfully.
%     IERR > 0  :  K steps were performed succesfully, but the algorithm
%                  switched to full reorthogonalization after IERR steps.
%     IERR < 0  :  Iteration was terminated after -IERR steps because an
%                  invariant subspace was found, and 3 deflation attempts
%                  were unsuccessful.
%   On exit WORK(1) contains the number of reorthogonalizations performed, and
%   WORK(2) contains the number of inner products performed in the
%   reorthogonalizations.
%

% References:
% R.M. Larsen, Ph.D. Thesis, Aarhus University, 1998.
%
% G. H. Golub & C. F. Van Loan, "Matrix Computations",
% 3. Ed., Johns Hopkins, 1996.  Chapter 9.
%
% B. N. Parlett, The Symmetric Eigenvalue Problem'',
% Prentice-Hall, Englewood Cliffs, NJ, 1980.
%
% H. D. Simon, The Lanczos algorithm with partial reorthogonalization'',
% Math. Comp. 42 (1984), no. 165, 115--142.

% Rasmus Munk Larsen, DAIMI, 1998

% Check input arguments.
if nargin<1, error('Not enough input arguments.');  end
if isnumeric(A) | isstruct(A)
if isnumeric(A)
[m n] = size(A);
if m~=n | ~isequal(A,A') | ~isreal(A)
error('A must be real symmetric')
end
elseif isstruct(A)
[m n] = size(A.L);
end

if nargin<7 | isempty(T_k),
anorm = []; est_anorm=1;
else
anorm = T_k; est_anorm=0;
end
if nargin<6,  Q_k=[]; T_k=[]; else,  T_k = Q_k; Q_k = options; end
if nargin<4 | isempty(r),  options = []; else,  options = r;  end
if nargin<3 | isempty(kmax),
r = rand(n,1)-0.5;
else
r = kmax;
end
if nargin<2 | isempty(nin);  kmax = max(10,n/10); else,  kmax = nin;  end
else
if nargin<2
error('Not enough input arguments.');
end
% Check input functions and parse to create an internal object
% if an explicit expression is given.
[A, msg] = fcnchk(A);
if ~isempty(msg)
error(msg);
end
n = nin;
if nargin<8 | isempty(anorm), anorm = []; est_anorm=1; else est_anorm=0; end
if nargin<7,  Q_k=[]; T_k=[]; end
if nargin<5 | isempty(options),  options = [];          end
if nargin<4 | isempty(r),  r = rand(n,1)-0.5;   end
if nargin<3 | isempty(kmax);  kmax = max(10,n/10); end
end

% Set options.
delta = sqrt(eps/kmax); % Desired level of orthogonality.
eta = eps^(3/4)/sqrt(kmax);     % Level of orth. after reorthogonalization.
cgs = 0;                % Flag for switching between iterated CGS and MGS.
elr = 1;                % Flag for switching extended local
% reorthogonalization on and off.
deflate = 0;              % Flag for deflation against OPTIONS.Y

% Parse options struct
if ~isempty(options) & isstruct(options)
c = fieldnames(options);
for i=1:length(c)
if strmatch(c(i),'delta'), delta = getfield(options,'delta');  end
if strmatch(c(i),'eta'), eta = getfield(options,'eta'); end
if strmatch(c(i),'cgs'), cgs = getfield(options,'cgs'); end
if strmatch(c(i),'elr'), elr = getfield(options,'elr'); end
if strmatch(c(i),'Y'), deflate = ~isempty(options.Y);  end
end
end

np = 0;  nr = 0; ierr=0;

% Rule-of-thumb estimate on the size of round-off terms:
eps1 = sqrt(n)*eps/2; % Notice that {\bf u} == eps/2.
gamma = 1/sqrt(2);

% Prepare Lanczos iteration
if isempty(Q_k) % New Lanczos tridiagonalization.
% Allocate space
alpha = zeros(kmax+1,1);  beta = zeros(kmax+1,1);
Q_k = zeros(n,kmax);
q = zeros(n,1); beta(1)=norm(r);
omega = zeros(kmax,1); omega_max = omega;  omega_old = omega;
omega(1) = 0;   force_reorth= 0;
j0 = 1;
else            % Extending existing Lanczos tridiagonalization.
j = size(Q_k,2); % Size of existing factorization
% Allocate space
Q_k = [Q_k zeros(n,kmax-j)];
alpha = zeros(kmax+1,1);  beta = zeros(kmax+1,1);
alpha(1:j) = diag(T_k);
if j>1
beta(2:j) = diag(T_k,-1);
end
q = Q_k(:,j);
% Reorthogonalize r.
beta(j+1) = norm(r);
if j<kmax & beta(j+1)*delta < anorm*eps1,
fro = 1;
end
if isfinite(delta)
int = 1:j;
[r,beta(j+1),rr] = reorth(Q_k,r,beta(j+1),int,gamma,cgs);
np = rr*j;    nr = 1;   force_reorth = 1;
else
force_reorth = 0;
end
% Compute Gerscgorin bound on ||T_k||_2 as SQRT(||T_k'*T_k||_1)
if est_anorm
anorm = sqrt(norm(T_k'*T_k,1));
end
omega = eps1*ones(kmax,1); omega_max = omega;  omega_old = omega;
j0 = j+1;
end

if delta==0
fro = 1; % The user has requested full reorthogonalization.
else
fro = 0;
end

for j=j0:kmax,
% Lanczos Step:
q_old = q;
if beta(j)==0
q = r;
else
q = r / beta(j);
end
Q_k(:,j) = q;
if isnumeric(A)
u = A*q;
elseif isstruct(A)
u = A.U \ ( A.L \ q);
else
u = feval(A,q);
end
r = u - beta(j)*q_old;
alpha(j) = q'*r;
r = r - alpha(j)*q;

% Extended local reorthogonalization:
beta(j+1) = sqrt(r'*r); % Quick and dirty estimate.
if beta(j+1)<gamma*beta(j) & elr
if  j==1
t1=0;
for i=1:2
t = q'*r;
r = r-q*t;
t1 = t1+t;
end
alpha(j) = alpha(j) + t1;
elseif j>1
t1 = q_old'*r;
t2 = q'*r;
r = r  - (q_old*t1 + q*t2); % Add small terms together first to
if beta(j)~=0               % reduce risk of cancellation.
beta(j) = beta(j) + t1;
end
alpha(j) = alpha(j) + t2;
end
beta(j+1) = sqrt(r'*r); % Quick and dirty estimate.
end

% Update Gersgorin estimate of ||T_k|| if required
%  if est_anorm & beta(j+1)~=0
%    T_k = spdiags([[beta(2:j);0] alpha(1:j) beta(1:j)],-1:1,j,j);
%    anorm = sqrt(norm(T_k'*T_k,1))
%  end
if  est_anorm & beta(j+1)~=0
anorm = update_gbound(anorm,alpha,beta,j);
end

% Update omega-recurrence
if j>1 & ~fro & beta(j+1)~=0
[omega,omega_old] = update_omega(omega,omega_old,j,alpha,beta,...
eps1,anorm);
omega_max(j) = max(abs(omega));
end

% Reorthogonalize if required
if j>1 & (fro  | force_reorth | omega_max(j)>delta) & beta(j+1)~=0
if fro
int = 1:j;
else
if force_reorth == 0
force_reorth= 1; % Do forced reorth to avoid spill-over from q_{j-1}.
int = compute_int(omega,j,delta,eta,0,0,0);
else
force_reorth= 0;
end
end
[r,beta(j+1),rr] = reorth(Q_k,r,beta(j+1),int,gamma,cgs);
omega(int) = eps1;
np = np + rr*length(int(:));    nr = nr + 1;
else
beta(j+1) = norm(r); % compute norm accurately.
end

if deflate
[r,beta(j+1),rr] = reorth(options.Y,r,beta(j+1),1:size(options.Y,2), ...
gamma,cgs);
end

if  j<kmax & beta(j+1) < n*anorm*eps  ,
% If beta is "small" we deflate by setting the off-diagonals of T_k
% to 0 and attempt to restart with a basis for a new
% invariant subspace by replacing r with a random starting vector:
beta(j+1) = 0;
bailout = 1;
for attempt=1:3
r = rand(n,1)-0.5;
if isnumeric(A)
r = A*r;
elseif isstruct(A)
r = A.U \ ( A.L \ r);
else
r = feval(A,r);
end
nrm=sqrt(r'*r); % not necessary to compute the norm accurately here.
int = 1:j;
[r,nrmnew,rr] = reorth(Q_k,r,nrm,int,gamma,cgs);
omega(int) = eps1;
np = np + rr*length(int(:));    nr = nr + 1;
if nrmnew > 0
% A vector numerically orthogonal to span(Q_k(:,1:j)) was found.
% Continue iteration.
bailout=0;
break;
end
end
if bailout
ierr = -j;
break;
else
r=r/nrmnew; % Continue with new normalized r as starting vector.
force_reorth = 1;
if delta>0
fro = 0;    % Turn off full reorthogonalization.
end
end
elseif j<kmax & ~fro & beta(j+1)*delta < anorm*eps1,
% If anorm*eps1/beta(j+1) > delta then  omega(j+1) will
% immediately exceed delta, and thus forcing a reorth. to occur at the
% next step. The components of omega will mainly be determined
% by the initial value and not the recurrence, and therefore we
% cannot tell reliably which components exceed eta => we might
% as well switch to full reorthogonalization to avoid trouble.
% The user is probably trying to determine pathologically
% small ( < sqrt(eps)*||A||_2 ) eigenvalues.
%    warning(['Semiorthogonality cannot be maintained at iteration ', ...
%	  num2str(j),'. The matrix is probably ill-conditioned.', ...
%	  ' Switching to full reorthogonalization.'])
fro = 1;
ierr = j;
end
end

% Set up tridiagonal T_k in sparse matrix data structure.
T_k = spdiags([[beta(2:j);0] alpha(1:j) beta(1:j)],-1:1,j,j);
if nargout<2
Q_k = T_k;
elseif j~=size(Q_k,2)
Q_k = Q_k(:,1:j);
end
work = [nr np];

function [omega,omega_old] = update_omega(omega, omega_old, j, ...
alpha,beta,eps1,anorm)
% UPDATE_OMEGA:  Update Simon's omega_recurrence for the Lanczos vectors.
%
% [omega,omega_old] = update_omega(omega, omega_old,j,eps1,alpha,beta,anorm)
%

% Rasmus Munk Larsen, DAIMI, 1998.

% Estimate of contribution to roundoff errors from A*v
%   fl(A*v) = A*v + f,
% where ||f|| \approx eps1*||A||.
% For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps.
T = eps1*anorm;
binv = 1/beta(j+1);

omega_old = omega;
% Update omega(1) using omega(0)==0.
omega_old(1)= beta(2)*omega(2)+ (alpha(1)-alpha(j))*omega(1) -  ...
beta(j)*omega_old(1);
omega_old(1) = binv*(omega_old(1) + sign(omega_old(1))*T);
% Update remaining components.
k=2:j-2;
omega_old(k) = beta(k+1).*omega(k+1) + (alpha(k)-alpha(j)).*omega(k) ...
+ beta(k).*omega(k-1) - beta(j)*omega_old(k);
omega_old(k) = binv*(omega_old(k) + sign(omega_old(k))*T);
omega_old(j-1) = binv*T;
% Swap omega and omega_old.
temp = omega;
omega = omega_old;
omega_old = omega;
omega(j) =  eps1;

function anorm = update_gbound(anorm,alpha,beta,j)
%UPDATE_GBOUND   Update Gerscgorin estimate of 2-norm
%  ANORM = UPDATE_GBOUND(ANORM,ALPHA,BETA,J) updates the Gerscgorin bound
%  for the tridiagonal in the Lanczos process after the J'th step.
%  Applies Gerscgorins circles to T_K'*T_k instead of T_k itself
%  since this gives a tighter bound.

if j==1 % Apply Gerscgorin circles to T_k'*T_k to estimate || A ||_2
i=j;
% scale to avoid overflow
scale = max(abs(alpha(i)),abs(beta(i+1)));
alpha(i) = alpha(i)/scale;
beta(i+1) = beta(i+1)/scale;
anorm = 1.01*scale*sqrt(alpha(i)^2+beta(i+1)^2 + abs(alpha(i)*beta(i+1)));
elseif j==2
i=1;
% scale to avoid overflow
scale = max(max(abs(alpha(1:2)),max(abs(beta(2:3)))));
alpha(1:2) = alpha(1:2)/scale;
beta(2:3) = beta(2:3)/scale;

anorm = max(anorm, scale*sqrt(alpha(i)^2+beta(i+1)^2 + ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2))));
i=2;
anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1))) );
elseif j==3
% scale to avoid overflow
scale = max(max(abs(alpha(1:3)),max(abs(beta(2:4)))));
alpha(1:3) = alpha(1:3)/scale;
beta(2:4) = beta(2:4)/scale;
i=2;
anorm = max(anorm,scale*sqrt(abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2))) );
i=3;
anorm = max(anorm,scale*sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1))) );
else
% scale to avoid overflow
%  scale = max(max(abs(alpha(j-2:j)),max(abs(beta(j-2:j+1)))));
%  alpha(j-2:j) = alpha(j-2:j)/scale;
%  beta(j-2:j+1) = beta(j-2:j+1)/scale;

% Avoid scaling, which is slow. At j>3 the estimate is usually quite good
% so just make sure that anorm is not made infinite by overflow.
i = j-1;
anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1) + alpha(i+1)*beta(i+1)) + ...
abs(beta(i+1)*beta(i+2)));
if isfinite(anorm1)
anorm = max(anorm,anorm1);
end
i = j;
anorm1 = sqrt(abs(beta(i)*beta(i-1)) + ...
abs(beta(i)*alpha(i-1) + alpha(i)*beta(i)) + ...
beta(i)^2+alpha(i)^2+beta(i+1)^2 +  ...
abs(alpha(i)*beta(i+1)));
if isfinite(anorm1)
anorm = max(anorm,anorm1);
end
end