function varargout = ahbeigs(varargin)
%
% AHBEIGS: will find a few eigenvalues and eigenvectors for either the
% standard eigenvalue problem A*x = lambda*x or the generalized eigenvalue
% problem A*x = lambda*B*x.
%
% [V,D,FLAG] = AHBEIGS(A,OPTIONS)
% [V,D,FLAG] = AHBEIGS(A,B,OPTIONS)
% [V,D,FLAG] = AHBEIGS('Afunc',N,B,OPTIONS)
%
% The first input argument must be a matrix A which can be passed as a numeric
% matrix or as a M-file ('Afunc') that computes the product A*X or inv(B)*A*X,
% where X is a (N x blsz) matrix. If A is passed as a M-file then the second
% input argument N is the size of the matrix A. For the generalized eigenvalue
% problem the matrix B is passed only as a numeric matrix. If a linear solver is
% required then use the M-file ('Afunc') and compute inv(B)*A*X. In all the
% implementations AHBEIGS(A,...) can be replaced with AHBEIGS('Afunc',N,...).
%
% The function 'Afunc' must be a M-file. The (N x blsz) matrix X is the first
% input argument, the second input argument is N, (the size of the matrix A),
% and the third input argument is blsz, i.e. Afunc(X,N,blsz). The output is
% the (N x blsz) matrix Y where, Y = A*X or Y = inv(B)*A*X.
%
% OUTPUT OPTIONS:
% ---------------
%
% I.) AHBEIGS(A) or AHBEIGS(A,B)
% Displays the desired eigenvalues.
%
% II.) D = AHBEIGS(A) or D = AHBEIGS(A,B)
% Returns the desired eigenvalues in the vector D.
%
% III.) [V,D] = AHBEIGS(A) or [V,D] = AHBEIGS(A,B)
% D is a diagonal matrix that contains the desired eigenvalues along the
% diagonal and the matrix V contains the corresponding eigenvectors, such
% that A*V = V*D or A*V = B*V*D. If AHBEIGS reaches the maximum number of
% iterations before convergence then the matrices V and D contain any
% eigenpairs that have converged plus the last Ritz pair approximations
% for the eigenpairs that have not converged.
%
% IV.) [V,D,FLAG] = AHBEIGS(A) or [V,D,FLAG] = AHBEIGS(A,B)
% This option returns the same as (III) plus a two dimensional array FLAG
% that reports if the algorithm converges and the number of matrix vector
% products. If FLAG(1) = 0 this implies normal return, all eigenvalues have
% converged. If FLAG(1) = 1 then the maximum number of iterations have been
% reached before all desired eigenvalues have converged. FLAG(2) contains the
% number of matrix vector products used by the code. If the maximum number of
% iterations are reached then the matrices V and D contain any eigenpairs
% that have converged plus the last Ritz pair approximations for the eigenpairs
% that have not converged.
%
% INPUT OPTIONS:
% --------------
%
% ... = AHBEIGS(A,OPTS) or ... = AHBEIGS(A,B,OPTS)
% OPTS is a structure containing input parameters. The input parameters can
% be given in any order. The structure OPTS may contain some or all of the
% following input parameters. The string for the input parameters can contain
% upper or lower case characters.
%
% INPUT PARAMETER DESCRIPTION
%
% OPTS.ADJUST Initial number of vectors to add to the K restart vectors. After
% vectors start to converge more vectors are added to help increase
% convergence. Should be set so that K + ADJUST is a multiple of BLSZ.
% DEFAULT VALUE OPTS.ADJUST = 3
%
% OPTS.BLSZ Block size of the Arnoldi Hessenberg matrix.
% DEFAULT VALUE BLSZ = 3
%
% OPTS.CHOLB Indicates if the Cholesky factorization of the matrix B is available.
% If the Cholesky factorization matrix R is available then set
% CHOLB = 1 and replace the input matrix B with R where B = R'*R.
% DEFAULT VALUE CHOLB = 0
%
% OPTS.DISPR Indicates if K Ritz values and residuals are to be displayed on each
% iteration. Set positive to display the Ritz values and residuals on
% each iteration.
% DEFAULT VALUE DISPR = 0
%
% OPTS.PERMB Permutation vector for the Cholesky factorization of B(PERMB,PERMB).
% When the input matrix B is replaced with R where B(PERMB,PERMB)=R'*R
% then the vector PERMB is the permutation vector.
% DEFAULT VALUE PERMB=1:N
%
% OPTS.K Number of desired eigenvalues.
% DEFAULT VALUE K = 6
%
% OPTS.MAXIT Maximum number of iterations, i.e. maximum number of restarts.
% DEFAULT VALUE MAXIT = 100
%
% OPTS.NBLS Number of blocks in the Arnoldi Hessenberg matrix. If number of
% blocks is not sufficiently large enough then AHBEIGS will not
% converge or miss some desired eigenvalues.
% DEFAULT VALUE NBLS = 10
%
% OPTS.SIGMA Two letter string or numeric value specifying the location
% of the desired eigenvalues.
% 'LM' or 'SM' Largest or Smallest Magnitude
% 'LR' or 'SR' Largest or Smallest Real part
% 'LI' or 'SI' Largest or Smallest Imaginary part
% 'LA' or 'SA' Largest or Smallest Algebraic (Symmetric Problems only)
% NVAL A numeric value. The program searches for the K closest
% eigenvalues to the numeric value NVAL. (Program will factor
% the matrix. This will increase storage requirements.)
% DEFAULT VALUE SIGMA = 'LM'
%
% OPTS.TOL Tolerance used for convergence. Convergence is determined when
% || Ax - lambda*x ||_2 <= TOL*max(abs(Ritz)).
% DEFAULT VALUE TOL = 1d-6
%
% OPTS.V0 A matrix of starting vectors.
% DEFAULT VALUE V0 = randn
%
% DATE: 6/15/2005
% VER: 1.0
% AUTHOR:
% James Baglama
% Department of Mathematics
% University of Rhode Island
% Kingston, Rhode Island 02881
% E-Mail: jbaglama@math.uri.edu
%
%
% Research supported by NSF grant DMS-0311786.
%
% Too many output arguments requested.
if (nargout >= 4) error('ERROR: Too many output arguments.'); end
%----------------------------%
% BEGIN: PARSE INPUT VALUES. %
%----------------------------%
% No input arguments, return help.
if nargin == 0, help ahbeigs, return, end
% Get matrix A. Check type (numeric or character) and dimensions.
if (isstruct(varargin{1})), error('A must be a matrix.'); end
A = varargin{1};
if ischar(varargin{1})
if nargin == 1, error('Need N (size of matrix A).'); end
n = varargin{2};
if ~isnumeric(n) | length(n) > 2
error('Second argument N must be a numeric value.');
end
else
[n,n] = size(varargin{1});
if any(size(varargin{1}) ~= n), error('Matrix A is not square.'); end
if ~isnumeric(varargin{1}), error('A must be a numeric matrix.'); end
end
% Set all input options to default values or set to empty arrays.
blsz = 3; adjust = []; cholB = 0; dispr = 0; K = 6; B=[];
maxit = 100; nbls = 10; permB = []; sigma = 'LM'; tol = 1d-6; V=[];
L=[]; U=[]; P=[]; nval=[];
% If a generalized eigenvalue problem is to be solved get matrix B.
% Check type (numeric or character) and dimensions.
if ((ischar(varargin{1}) & nargin > 2) | (isnumeric(varargin{1}) & nargin > 1))
if ~isstruct(varargin{2+ischar(varargin{1})})
B = varargin{2+ischar(varargin{1})};
if ~isnumeric(B), error('B must be a numeric matrix.'); end
if any(size(B) ~= n), error('Matrix B must be the same size as A.'); end
end
end
% Get input options from the structure array.
if nargin > 1 + ischar(varargin{1}) + ~isempty(B)
options = varargin{2+ischar(varargin{1}) + ~isempty(B):nargin};
names = fieldnames(options);
I = strmatch('ADJUST',upper(names),'exact');
if ~isempty(I), adjust = getfield(options,names{I}); end
I = strmatch('BLSZ',upper(names),'exact');
if ~isempty(I), blsz = getfield(options,names{I}); end
I = strmatch('CHOLB',upper(names),'exact');
if ~isempty(I), cholB = getfield(options,names{I}); end
I = strmatch('DISPR',upper(names),'exact');
if ~isempty(I), dispr = getfield(options,names{I}); end
I = strmatch('K',upper(names),'exact');
if ~isempty(I), K = getfield(options,names{I}); end
I = strmatch('MAXIT',upper(names),'exact');
if ~isempty(I), maxit = getfield(options,names{I}); end
I = strmatch('NBLS',upper(names),'exact');
if ~isempty(I), nbls = getfield(options,names{I}); end
I = strmatch('PERMB',upper(names),'exact');
if ~isempty(I), permB = getfield(options,names{I}); end
I = strmatch('SIGMA',upper(names),'exact');
if ~isempty(I), sigma = upper(getfield(options,names{I})); end
I = strmatch('TOL',upper(names),'exact');
if ~isempty(I), tol = getfield(options,names{I}); end
I = strmatch('V0',upper(names),'exact');
if ~isempty(I), V = getfield(options,names{I}); end
end
% Check type of input values and output error message if needed.
if (~isnumeric(blsz) | ~isnumeric(K) | ~isnumeric(dispr) | ...
~isnumeric(maxit) | ~isnumeric(nbls) | ~isnumeric(tol) | ...
~isnumeric(permB) | ~isnumeric(cholB))
error('Incorrect type for input value(s) in the structure.');
end
% If a numeric value is given for sigma then set nval=sigma and
% sigma = 'LM' to denote that the code is searching eigenvalues near nval.
% Using the transformation inv(A-sigma*B)*B x = theta * x where
% theta = 1/(lambda-sigma).
if isnumeric(sigma), nval = sigma; sigma = 'LM'; end
% Check the length and value of the character sigma.
if length(sigma) ~= 2, error('SIGMA must be SA, LA, SM, LM, SR, LR, LI, or SI.'), end
if (~strcmp(sigma,'SA') & ~strcmp(sigma,'LA') & ~strcmp(sigma,'LI') & ...
~strcmp(sigma,'SM') & ~strcmp(sigma,'LM') & ~strcmp(sigma,'SI') & ...
~strcmp(sigma,'SR') & ~strcmp(sigma,'LR'))
error ('SIGMA must be SA, LA, SM, LM, SR, LR LI, or SI.');
end
% Resize Krylov subspace if blsz*nbls (i.e. number of Arnoldi vectors) is larger
% than n (i.e. the size of the matrix A).
if blsz*nbls >= n
nbls = floor(n/blsz);
warning(['Changing NBLS to ',num2str(nbls)]);
end
% Increase the number of desired values to help increase convergence.
% Set K+adjust to be next multiple of BLSZ > K.
if isempty(adjust)
adjust = 1:blsz; I = find(mod(K + adjust-1,blsz) == 0); adjust = adjust(I(1));
end
K_org = K; K = K + adjust; % K_org to be orginal value of K.
% Check for input errors in the structure array.
if K <= 0, error('K must be a positive value.'); end
if K > n, error('K is too large. Use eig(full(A))! '); end
if blsz <= 0, error('BLSZ must be a positive value.'); end
if nbls <=1, error('NBLS must be greater than 1.'); end
if tol < 0, error('TOL must be non-negative.'); end
if maxit <=0, error('MAXIT must be positive.'); end
if ((cholB ~= 0) & (cholB ~= 1)), error('Unknown value for CHOLB.'); end
if ~isempty(permB)
if ((size(permB,1) ~= n) | (size(permB,2) ~= 1)) & ...
((size(permB,1) ~= 1) | (size(permB,2) ~= n))
error('Incorrect size for PERMB.');
end
end
% Automatically adjust Krylov subspace to accommodate larger values of K.
if blsz*(nbls-1) - blsz - K - 1 < 0
nbls = ceil((K+1)/blsz+2.1);
warning(['Changing NBLS to ',num2str(nbls)]);
end
if blsz*nbls >= n,
nbls = floor(n/blsz-0.1);
warning(['Changing NBLS to ',num2str(nbls)]);
end
if blsz*(nbls-1) - blsz - K - 1 < 0, error('K is too large.'); end
% If starting matrix V0 is not given then set starting matrix V0 to be a
% (n x blsz) matrix of normally distributed random numbers. Also,
% preallocate memory size, so no resizing of the large matrix V occurs.
if nnz(V) == 0
V = [randn(n,blsz) zeros(n,nbls*blsz)];
else
if ~isnumeric(V), error('Incorrect starting matrix V0.'); end
if (size(V,1) ~= n), error('Incorrect size of starting matrix V0.'); end
if nnz(V) == 0, error('Matrix V contains all zeros.'); end
V = [V zeros(n,nbls*blsz)];
end
% Set tolerance to machine precision if tol < eps.
if tol < eps, tol = eps; end
% Using shift and invert if sigma is given as a real number.
if ~isempty(nval)
if ~isnumeric(varargin{1})
error('A must be a numeric matrix in order to use the internal shift and invert.');
end
if ~isempty(B)
[L,U,P] = lu(varargin{1}-nval*B);
else
[L,U,P] = lu(varargin{1}-nval*eye(n));
end
end
% If a generalized eigenvalue problem is to be solved then try to get the Cholesky
% factorization of the matrix B (i.e. B=R'*R) and set B=R the upper triangular matrix.
if ~isempty(B) & isempty(nval)
% If B is sparse then compute the sparse cholesky factorization. Using symmmd or
% symamd to get a permutation vector permB such that B(permB,permB) will have a
% sparser cholesky factorization than B. The permutation is given by
% B(permB,permB) = I(:,permB)'*B*I(:,permB), where I is the identity matrix.
% To compute the product A(permB,permB)*X use Y(permB,:) = X; Y=A*Y; Y=Y(permB,:);
if isempty(permB), permB = 1:n; end
% If cholB = 1, then the input matrix B is given as R in the Cholesky factorization.
if ~cholB
if ~isreal(B), error('Matrix B must be real.'); end
[L,cholerr] = chol(B(permB,permB));
% Will try to solve problem using the inverse of B, (i.e. inv(B)*A*x = lambda*x)
% by computing the LU factorization of B. This may cause numerical instabilities.
if cholerr
[L,U,P] = lu(B); % Used only if B is not symmetric positive definite.
else
B = L; L=[];
end
end
end
%--------------------------%
% END: PARSE INPUT VALUES. %
%--------------------------%
%-----------------------------------------------------------%
% BEGIN: DESCRIPTION AND INITIALIZATION OF LOCAL VARIABLES. %
%-----------------------------------------------------------%
check_norm_balance=[]; % Accuray of eigenpairs of H. Determines if nobalance in eig is needed.
conv = 'F'; % Boolean to determine if all desired eigenvalues have converged.
H = []; % Block Hessenberg (Hsz_n x Hsz_m) matrix.
H_R = []; % Stores the last block of H. Used in updating H.
H_eig = []; % Eigenvalues of H.
H_eig_no = []; % Eigenvalues of H using nobalance in eig. Not always computed.
H_eigv = []; % Eigenvectors of H.
H_eigv_no = []; % Eigenvectors of H using nobalance in eig. Not always computed.
Hsz_m = []; % Second dimension of H.
Hsz_n = []; % First dimension of H.
iter = 1; % Main loop iteration count.
m = blsz; % Current number of columns in the matrix V.
mprod = 0; % The number of matrix vector products.
Q = []; % Schur vectors for H.
R = []; % Is a diagonal matrix of 1 or -1's used in the function setuparnv.
sqrteps = sqrt(eps); % Square root of machine tolerance.
T = []; % T matrix in the WY representation of the Householder products.
T_schur=[]; % Schur matrix for H.
Tsz_m = []; % Size of T.
%---------------------------------------------------------%
% END: DESCRIPTION AND INITIALIZATION OF LOCAL VARIABLES. %
%---------------------------------------------------------%
%-----------------------------%
% BEGIN: MAIN ITERATION LOOP. %
%-----------------------------%
while (iter <= maxit)
% Compute the block Householder Arnoldi decomposition.
[V,H,T,mprod] = ArnHouBlk(varargin{1},B,L,U,P,nval,permB,V,H,T,n,m,nbls,blsz,mprod);
% Determine the size of the block Hessenberg matrix H. Possible
% truncation may occur in ArnHouBlk if an invariant subspace
% has been found or K + ADJUST is not a multiple of BLSZ.
Hsz_n = size(H,1); Hsz_m = size(H,2);
% Compute the eigenvalue decomposition of the block Hessenberg H(1:Hsz_m,:).
[H_eigv,H_eig] = eig(H(1:Hsz_m,:));
% Check the accuracy of the computation of the eigenvalues of the
% Hessenberg matrix. This is used to monitor balancing.
check_norm_balance = norm(H(1:Hsz_m,:)*H_eigv - H_eigv*H_eig);
if check_norm_balance >= eps*norm(H(1:Hsz_m,:))
[H_eigv_no,H_eig_no] = eig(H(1:Hsz_m,:),'nobalance');
check_norm_nobalance = norm(H(1:Hsz_m,:)*H_eigv_no - H_eigv_no*H_eig_no);
if check_norm_nobalance < check_norm_balance
H_eigv = H_eigv_no; H_eig = H_eig_no;
end
end
H_eig = diag(H_eig);
% Sort the eigenvalues and check for convergence.
[H_eig,H_eigv,K,conv] = eigsortconv(Hsz_m,Hsz_n,blsz,dispr,H,iter,K_org,...
K,adjust,H_eig,H_eigv,sigma,tol);
% If all desired Ritz values converged then exit main loop.
if strcmp(conv,'T'), break; end
% Update the main iteration loop count.
iter = iter + 1; if iter >= maxit, break; end
% Convert the complex eigenvectors of the eigenvalue decomposition of H
% to real vectors and convert the complex diagonal matrix to block diagonal.
% Use MATLAB's internal function cdf2rdf to accomplish this task.
if ~isreal(H_eig)
[Q,T_Schur] = cdf2rdf(H_eigv,diag(H_eig)); % Convert complex vectors to real vectors.
[Q,R] = qr(Q(:,1:K),0); % Compute the QR factorization of H_eigv(:,1:K).
else
[Q,R] = qr(H_eigv(:,1:K),0); % Compute the QR factorization of H_eigv(:,1:K).
end
% The Schur matrix for H.
T_Schur = triu(Q'*H(1:Hsz_m,:)*Q,-1);
% Compute the starting vectors and the residual vectors from the Householder
% WY form. The starting vectors will be the first K Schur vectors and the
% residual vectors are stored as the last BLSZ vectors in the Householder WY form.
Tsz_m = size(T,1);
V(:,Hsz_n-blsz+1:Hsz_n)= V(:,1:Tsz_m)*(T*V(Tsz_m-blsz+1:Tsz_m,1:Tsz_m)');
V(Tsz_m-blsz+1:Tsz_m,Hsz_n-blsz+1:Hsz_n) = eye(blsz,blsz)+V(Tsz_m-blsz+1:Tsz_m,Hsz_n-blsz+1:Hsz_n);
%
V(:,1:K) = V(:,1:Hsz_m)*(T(1:Hsz_m,1:Hsz_m)*(V(1:Hsz_m,1:Hsz_m)'*Q(:,1:K)));
V(1:Hsz_m,1:K) = Q(:,1:K) + V(1:Hsz_m,1:K);
% Set the size of the large matrix V and move the residual vectors.
m = K + 2*blsz; V(:,K+1:K+blsz) = V(:,Hsz_n-blsz+1:Hsz_n);
% Set the new starting vector(s) to be the desired vectors V(:,1:K) with the
% residual vectors V(:,Hsz_n-blsz+1:Hsz_n). Place all vectors in the compact
% WY form of the householder product. Compute the next set of vectors by
% computing A*V(:,Hsz_n-blsz+1:Hsz_n) and store this in V(:,Hsz_n+1:Hsz_n+blsz).
[V,T,R,mprod] = SetupArnVecs(V,T,varargin{1},B,L,U,P,nval,permB,n,blsz,m-blsz,mprod);
% Compute the first K columns and K+BLSZ rows of the matrix H, used in augmenting.
H_R = H(Hsz_n-blsz+1:Hsz_n,Hsz_m-blsz+1:Hsz_m); H=zeros(K+blsz,K);
H(1:K,1:K) = R(1:K,1:K)*T_Schur(1:K,1:K)*R(1:K,1:K);
H(K+1:K+blsz,1:K) = R(K+1:K+blsz,K+1:K+blsz)*H_R*Q(Hsz_m-(blsz-1):Hsz_m,1:K)*R(1:K,1:K);
end
%---------------------------%
% END: MAIN ITERATION LOOP. %
%---------------------------%
%------------------------%
% BEGIN: OUTPUT RESULTS. %
%------------------------%
% Test to see if maximum number of iterations have been reached.
FLAG(1) = 0; if iter >= maxit FLAG(1) = 1; end;
if FLAG(1)==1, warning('AHBEIGS did not find all requested eigenvalues!'); end
% Truncated eigenvalue and eigenvector arrays to include only desired eigenpairs.
H_eig = H_eig(1:K_org); H_eigv = H_eigv(:,1:K_org);
% Convert the shift and invert eigenvalues back to orginal state using 1/(lambda-sigma).
if ~isempty(nval), H_eig = 1./H_eig + nval; mprod = []; end
% Sort output.
[sortval,J] = sort(abs(H_eig(1:K_org))); J = J(length(J):-1:1);
% Output option I: Display eigenvalues only.
if (nargout == 0), eigenvalues = H_eig(J), end
% Output option II: Set eigenvalues equal to output vector.
if (nargout == 1), varargout{1} = H_eig(J); end
% Output option III. Output diagonal matrix of eigenvalues and
% corresponding matrix of eigenvectors.
if nargout == 2 | nargout == 3
H_eigv = H_eigv(:,1:K_org); H_eigv = H_eigv(:,J);
Tsz_m = size(T,1); H_sz1 = size(H_eigv,1); H_sz2 = size(H_eigv,2);
V = V(:,1:Tsz_m)*(T*(V(1:H_sz1,1:Tsz_m)'*H_eigv));
V(1:H_sz1,1:H_sz2) = H_eigv + V(1:H_sz1,1:H_sz2);
% Must solve a linear system to extract generalized eigenvectors.
if ~isempty(B) & isempty(L) & isempty(nval), V = B\V(permB,:); end
varargout{1} = V; varargout{2} = diag(H_eig(J));
end
% Output option IV: Output FLAG.
if nargout == 3, FLAG(2)=mprod; varargout{3}=FLAG; end
matrix_prod = mprod
%----------------------%
% END: OUTPUT RESULTS. %
%----------------------%
%---------------------------------------------%
% BEGIN: AUGMENTED BLOCK HOUSEHOLDER ARNOLDI. %
%---------------------------------------------%
function [V,H,T,mprod] = ArnHouBlk(A,B,L,U,P,nval,permB,V,H,T,n,m,nbls,blsz,mprod)
% Computes the Block Arnoldi decomposition using Householder reflections.
%
% The matrix A can be passed as a numeric matrix or as a filename. Note
% that if the matrix A is a filename then the file must accept [X,N,BLSZ]
% as input in that order and the file must return the matrix-vector product A*X.
%
% Input:
% A - Matrix A.
% B - Matrix B for generalized eigenvalue problems.
% L - LU factorization of (A-nval*B) or B.
% U - LU factorization of (A-nval*B) or B.
% P - Permutation matrix for the LU factorization.
% NVAL - Numeric value of sigma, if given.
% PERMB - Permutation vector for generalized eigenvalue problems.
% V - (N x M) Initial matrix to start the Arnoldi reduction.
% H - Empty or contains Schur values and residuals values on restart.
% T - Upper triangular matrix used in Householder WY representation.
% N - Size of the matrix A.
% M - Current size of V.
% NBLS - Maximum number of Arnoldi vectors.
% BLSZ - Block size.
% MPROD - Number matrix-vector products.
%
% Output:
% V - N x BLSZ*NBLS+BLSZ matrix in WY-Householder format.
% H - NBLS*BLSZ+BLSZ x NBLS*BLSZ upper block Hessenberg matrix.
% T - Upper triangular matrix used in Householder WY representation.
% MPROD - Number of matrix-vector products.
% James Baglama
% DATE: 7/20/2005
% LOCAL VARIABLES
T_blk =zeros(blsz,blsz); % Small BLSZ x BLSZ upper triangular matrix used in Householder WY representation.
% Check size of the block Arnoldi matrix H. If size is greater than n
% reduce size.
if blsz*(nbls+1) > n, nbls = floor(n/blsz)-1; end
% Begin of main iteration loop for the Augmented Block Householder Arnoldi decomposition.
while (m <= blsz*(nbls+1))
% Compute Householder vectors for BLSZ vectors V(m-blsz+1:n,m-blsz+1:m) using
% MATLAB's internal qr algorithm.
V(m-blsz+1:n,m-blsz+1:m) = qr(V(m-blsz+1:n,m-blsz+1:m));
% Compute the columns of the upper block hessenberg matrix H.
if m > blsz
H(1:m-blsz,m-2*blsz+1:m-blsz) = V(1:m-blsz,m-blsz+1:m);
H(m-blsz+1:m,m-2*blsz+1:m-blsz) = triu(V(m-blsz+1:m,m-blsz+1:m));
end
% Set up householder vectors.
m2 = 2*m; if m2 > n, m2=n; end
V(1:m2,m-blsz+1:m) = tril(V(1:m2,m-blsz+1:m),-(m-blsz+1));
for i=m-blsz+1:m,V(i,i)=1; end
% Using the compact WY form of the householder product. See,
% "A storage-efficient WY representation for products of
% householder transformations", by Schreiber and Van Loan,
% SIAM J. Sci. Stat. Comput. Vol. 10, No. 1 (1989) pp. 53-57.
% P_{1} * P_{2} * ... * P_{n} = I + V *T *V';
% Place BLSZ householder vectors into I + V *T_blk *V' form.
T_blk(1,1) = -2/(V(:,m-blsz+1)'*V(:,m-blsz+1));
for i=2:blsz
T_blk(1:i,i) = (V(:,m-blsz+i)'*V(:,m-blsz+1:m-blsz+i))';
T_blk(i,i) = -2/T_blk(i,i);
T_blk(1:i-1,i) = T_blk(i,i)*T_blk(1:i-1,1:i-1)*T_blk(1:i-1,i);
end
% Update I + V *T *V'.
T = [ [T; zeros(blsz,m-blsz)] [(T*(V(:,m-blsz+1:m)'*V(:,1:m-blsz))')*T_blk; T_blk]];
if m <= blsz*nbls
% V_{m+blsz} = ( I + V_m * T_m * V_m') * E_m.
V(:,m+1:m+blsz) = V(:,1:m)*(T*V(m-blsz+1:m,1:m)');
V(m-blsz+1:m,m+1:m+blsz) = eye(blsz,blsz)+V(m-blsz+1:m,m+1:m+blsz);
% Matrix-vector product, V_{m+blsz} = A * V_{m+blsz}.
[V(:,m+1:m+blsz),mprod] = matrixprod(A,B,L,U,P,nval,permB,V(:,m+1:m+blsz),n,blsz,mprod);
% V_{m+blsz} = (I + V_m * T_m' * V_m')*V_{m+blsz}.
V(:,m+1:m+blsz) = V(:,m+1:m+blsz) + V(:,1:m)*((V(:,m+1:m+blsz)'*V(:,1:m))*T)';
end
% Update iteration count.
m = m + blsz;
end % End main loop.
%-------------------------------------------%
% END: AUGMENTED BLOCK HOUSEHOLDER ARNOLDI. %
%-------------------------------------------%
%--------------------------------------------------------%
% BEGIN: EIGENVALUE REORDERING AND CONVERGENCE TESTING. %
%--------------------------------------------------------%
function [H_eig,H_eigv,K,conv] = eigsortconv(Hsz_m,Hsz_n,blsz,dispr,H,iter,K_org,K,...
adjust,H_eig,H_eigv,sigma,tol);
% This reorders the eigenvalues of H and determines convergence.
%
% Input:
% HSZ_M - Second dimension of H.
% HSZ_N - First dimension of H.
% BLSZ - Block size of H.
% DISPR - Boolean to determine if Ritz pairs should be displayed.
% H - Block Hessenberg (Hsz_n x Hsz_m) matrix.
% ITER - Number of iterations of the main loop.
% K_ORG - Orginal number of desired eigenpairs.
% K - Adjusted number of eigenpairs.
% ADJUST - Orginal value added to K_ORG
% H_EIG - Eigenvalues of H.
% H_EIGV - Eigenvectors of H.
% SIGMA - Two letter string specifying the location of the desired eigenvalues.
% TOL - Tolerance used for convergence. Convergence is determined when
% || Ax - lambda_i*x ||_2 <= TOL*abs(Ritz_i).
%
% Output:
% H_EIG - Sorted eigenvalues of H.
% H_EIGV - Sorted eigenvectors of H.
% K - Adjusted number of eigenpairs. Increased if eigenvectors converge.
% CONV - Boolean to indicate all eigenpairs have converged.
% James Baglama
% DATE: 7/20/2005
% LOCAL VARIABLES
conv = 'F'; % Boolean to determine if all desired eigenpairs have converged.
eig_conj=[]; % Indicator for when K splits a conjugate pair. K is increased by 1.
I_eig = []; % Integer vector of indexes for the ordering of the eigenvalues.
Len_res=[]; % Number of converged values.
% If sigma = LI or SI and all eigenvalues are real then sort by the real
% value.
if isreal(H_eig) & (strcmp(sigma,'SI') | strcmp(sigma,'LI'))
if strcmp(sigma,'SI'), sigma = 'SR'; end
if strcmp(sigma,'LI'), sigma = 'LR'; end
end
% Determine the ordering of the eigenvalues.
switch sigma
% Eigenvalues of largest magnitude.
case 'LM', [Y,I_eig] = sort(abs(H_eig)); I_eig = I_eig(length(I_eig):-1:1);
% Eigenvalues of smallest magnitude.
case 'SM', [Y,I_eig] = sort(abs(H_eig));
% Eigenvalues of largest real part.
case 'LR', [Y,I_eig] = sort(real(H_eig)); I_eig = I_eig(length(I_eig):-1:1);
% Eigenvalues of smallest real part.
case 'SR', [Y,I_eig] = sort(real(H_eig));
% Eigenvalues of largest imaginary part.
case 'LI', [Y,I_eig] = sort((abs(imag(H_eig)))); I_eig = I_eig(length(I_eig):-1:1);
% Eigenvalues of smallest imaginary part.
case 'SI', [Y,I_eig] = sort((abs(imag(H_eig))));
% Largest algebraic eigenvalues (symmetric problems only).
case 'LA', [Y,I_eig] = sort(H_eig); I_eig = I_eig(length(I_eig):-1:1);
% Smallest algebraic eigenvalues (symmetric problems only).
case 'SA', [Y,I_eig] = sort(H_eig);
end
% Sort the eigenvalue and eigenvector arrays.
H_eig = H_eig(I_eig); H_eigv = H_eigv(:,I_eig);
% Compute the residuals for the K_org Ritz values.
residuals = (sum(abs(H(Hsz_n-blsz+1:Hsz_n,Hsz_m-blsz+1:Hsz_m)*...
H_eigv(Hsz_m-(blsz-1):Hsz_m,1:K_org)).^2, 1).^(0.5));
% Output intermediate results.
if dispr ~= 0
format short e
% Sort output values.
[sortval,J] = sort(abs(H_eig(1:K_org))); J = J(length(J):-1:1);
S = [H_eig(J) residuals(J)'];
if ~isreal(S)
disp(sprintf(' Ritz values Residual Iteration: %d',iter));
else
disp(sprintf(' Ritz values Residual Iteration: %d',iter));
end
disp(S); disp(' '); disp(' ');
end
% Check for convergence.
Len_res = length(find(residuals(1:K_org)' < tol*abs(H_eig(1:K_org))));
if Len_res >= K_org, conv = 'T'; end % Set convergence to true.
% Adjust K to include more vectors as the number of vectors converge.
Len_res = length(find(residuals(1:K_org)' < eps*abs(H_eig(1:K_org))));
K = K_org + adjust+ Len_res; if K > Hsz_m - 2*blsz-1, K = Hsz_m - 2*blsz-1; end
% Determine if K splits a conjugate pair. If so replace K with K + 1.
if ~isreal(H_eig(K))
eig_conj = 1;
if K < Hsz_m
if abs(imag(H_eig(K)) + imag(H_eig(K+1))) < sqrt(eps)
K = K + 1; eig_conj = 0;
end
end
if K > 1 & eig_conj
if abs(imag(H_eig(K)) + imag(H_eig(K-1))) < sqrt(eps)
eig_conj = 0;
end
end
if eig_conj, error('%d th conjugate pair split.',K(i)); end
end
%-----------------------------------------------------%
% END: EIGENVALUE REORDERING AND CONVERGENCE TESTING. %
%-----------------------------------------------------%
%----------------------------------------------------------%
% BEGIN: WY HOUSEHOLDER REPRESENATION OF STARTING VECTORS. %
%----------------------------------------------------------%
function [V,T,R,mprod] = SetupArnVecs(V,T,A,B,L,U,P,nval,permB,n,blsz,m,mprod)
% Given an orthogonal matrix V, this function SetupArnVecs computes the
% compact WY form of the householder product. See, "A storage-efficient
% WY representation for products of householder transformations", by Schreiber
% and Van Loan, SIAM J. Sci. Stat. Comput. Vol. 10, No. 1 (1989) pp. 53-57.
% Also to save on computations the next vector in the block Arnoldi process
% is also calculated.
%
% Input:
% V - (N x M) orthogonal Matrix.
% T - Upper triangular matrix for WY-compact form.
% A - Matrix A.
% B - Matrix B for generalized eigenvalue problems.
% L - LU factorization of (A-nval*B) or B.
% U - LU factorization of (A-nval*B) or B.
% P - Permutation matrix for the LU factorization.
% NVAL - Numeric value of sigma, if given.
% PERMB - Permutation vector for generalized eigenvalue problems.
% N - Size of the matrix A.
% BLSZ - Block size
% m - Current size of V.
% MPROD - Number of matrix-vector products.
%
% Output:
% V - (N x M+BLSZ) matrix such that (N x 1:M) are the Householder vectors
% and (N x M+1:M+BLSZ) are the Arnoldi vectors needed to restart the
% Arnoldi process
% T - (M x M) Upper triangular matrix such that P_{1} * ... * P_{n} = I + V*T *V'.
% R - R is a diagonal matrix of 1's and -1's. R*R = I where
% V(input) * R = (I + V * T * V') * I(:,m).
% MPROD - Number of matrix-vector products.
% James Baglama
% DATE: 7/20/2005
% LOCAL VARIABLES
V_sign = []; % Used to determine Householder vectors, the sign of the diagonal value.
Vdot = []; % Dot product of Householder vectors.
D = []; % Used for scaling.
% Initialize matrix T.
T = V(1:m,1:m);
% Matrix-vector product, V_{m+1:m+blsz} = A * V_{m-blsz+1:m}.
[V(:,m+1:m+blsz),mprod] = matrixprod(A,B,L,U,P,nval,permB,V(:,m-blsz+1:m),n,blsz,mprod);
% Compute the householder vectors V. The input matrix V(:,1:m) is orthogonal
% and the length of each vector is 1. Hence the householder vector
% u_i = v_i + e_i*||v_i||_2 *sign(v_ii) = v_i + e_i*sign(v_ii) and the
% dot product u_i'*u_i = 2 + 2*sign(v_ii)*v_ii. Also, since the matrix is
% orthognal only the first m columns and rows are updated at this stage.
% At this stage V=[V_1; V_2], V_1 is a m x m matrix which is now computed.
for i =1:m
V_sign(i) = sign(V(i,i)); if V_sign(i) == 0, V_sign(i)=1; end
Vdot = 1 + V_sign(i)*V(i,i); % Dot product of householder vectors.
V(i,i) = V(i,i) + V_sign(i); % Reflection to the ith axis.
D(i,i) = 1/V(i,i); % Used for scaling. Note: V(i,i) >= 1.
% Apply the householder vectors to update the i+1:m vectors in V. Only
% the first m rows are needed at this stage.
V(1:m,i+1:m) = V(1:m,i+1:m) - V(1:m,i)*((V_sign(i)/Vdot)*V(i,i+1:m));
end
% Scale the matrix V_1, so that diagonal elements are 1.
V(1:m,1:m) = V(1:m,1:m)*D;
% Compute the upper triangular matrix R such that the following relationship
% holds V(input) * R = (I + V * T * V') * I(:,m).
R = -(diag(V_sign));
% Update the new block.
V(:,m+1:m+blsz) = V(:,m+1:m+blsz)*(R(m-blsz+1:m,m-blsz+1:m));
% Set up householder vectors. Put exact zeros in place of computed zeros.
V(1:m,1:m) = tril(V(1:m,1:m),0);
% Compute T. Note that V(input) * R = (I + V * T * V') * E_m. Let
% V = [V_1; V_2]. V_1 is already computed. Hence we have
% T = inv(V_1)*(E_m'*V_(input)*R - E_m)*inv(V_1)'.
% V_1 is upper triangular matrix with 1's on the diagonal.
T = T*R - eye(m); T = T/V(1:m,1:m)'; T = triu(V(1:m,1:m)\T);
% Compute the V_2 n-m x m matrix in V = [V_1; V_2].
% V_2 = V(input)*[0; I(n-m,m)]*R*inv(T*V_1).
V(m+1:n,1:m) = V(m+1:n,1:m)*(R/(T*V(1:m,1:m)'));
% V_{m+blsz} = (I + V_m * T_m' * V_m')*V_{m+blsz}. Update the vectors to
% be used in the block Arnoldi subroutine.
V(:,m+1:m+blsz) = V(:,m+1:m+blsz) + V(:,1:m)*((V(:,m+1:m+blsz)'*V(:,1:m))*T)';
%---------------------------------------------------------%
% END: WY HOUSEHOLDER REPRESENATION OF STARTING VECTORS. %
%---------------------------------------------------------%
%-------------------------------%
% BEGIN: MATRIX-VECTOR PRODUCT. %
%-------------------------------%
function [X,mprod] = matrixprod(A,B,L,U,P,nval,permB,X,n,blsz,mprod)
% Computes the matrix vector products.
%
% Input:
% A - Matrix A.
% B - Matrix B for generalized eigenvalue problems.
% L - LU factorization of (A-nval*B) or B.
% U - LU factorization of (A-nval*B) or B.
% P - Permutation matrix for the LU factorization.
% NVAL - Numeric value of sigma, if given.
% PERMB - Permutation vector for generalized eigenvalue problems.
% X - (N x BLSZ) Matrix to multiply OP (operator).
% N - Size of the matrix A.
% BLSZ - Block size.
% MPROD - Number of matrix-vector products.
%
% Output:
% X - (N x BLSZ) Product of OP*X (operator).
% MPROD - Number of matrix-vector products.
% James Baglama
% DATE: 7/20/2005
% Generalized eigenvalue problem (Cholesky factored).
if ~isempty(B) & isempty(L) & isempty(nval)
X(permB,:) = B \ X;
if ischar(A), X = feval(A,X,n,blsz); else X = A*X; end
X = (X(permB,:)'/B)';
end
% Generalized eigenvalue problem (LU factored).
if ~isempty(B) & ~isempty(L) & isempty(nval)
if ischar(A), X = feval(A,X,n,blsz); else X = A*X; end
X = U \ (L \ (P * X));
end
% Generalized or standard eigenvalue problem (Shift & Invert).
if ~isempty(nval)
if isempty(B)
X = U \ (L \ (P * X));
else
X = U \ (L \ (P * B* X));
end
end
% Standard eigenvalue problem.
if isempty(B) & isempty(nval)
if ischar(A), X = feval(A,X,n,blsz); else X = A*X; end
end
% Number of matrix-vector products with A.
mprod = mprod+blsz;
%-----------------------------%
% END: MATRIX-VECTOR PRODUCT. %
%-----------------------------%