function ShowEigs()
randn('seed',0)
rand('seed',0)
% Illustrates various ways to use the Matlab sparse eigenvalue solver eigs.
%
% We use the n-by-n arrow matrix
%
% d(1) 0 0 0 v(1)
% 0 d(2) 0 0 v(2)
% A = 0 0 d(3) 0 v(3) (n=5)
% 0 0 0 d(4) v(4)
% w(1) w(2) w(3) w(4) d(5)
%
% as an example.

% 1. The main user-responsibility when using eigs is to figure out how
% to organize matrix-vector products. One way is to pass eigs the matrix
% in sparse format. This requires a matrix set-up...

d = rand(4,1); v = rand(3,1); w = rand(3,1); A = sparseSetUp(d,v,w);

% 2. Alternatively, you can pass to eigs a function that carries out the
% operation vector = A*vector. This is preferable if the sparsity is
% structured and nicely parameterized.

x = rand(4,1); y = MyProd(x,d,v,w); yAnotherWay = A*x;

% Eigs can be called in several ways. We illustrate some typical cases
% below. Type "help eigs" for more detail.

% Data...
n = 20; d = (n:-1:1)'; v = n*randn(n-1,1); w = n*randn(n-1,1);
A = sparseSetUp(d,v,w);

% 3. Two ways to get the m largest in magnitude ('LM') eigenvalues and eigenvectors
% with default tolerances...
m = 3;
[V,D] = eigs(A,m,'LM');
[V,D] = eigs(@(x) MyProd(x,d,v,w),n,m,'LM');
% A*V(:,k) = D(k,k)*V(:,k) (approximately) for k=1:m.

% 4. Other types of extremal eigenvalues are possible. Eg., smallest in magnitude,
% largest or smallest real part, largest or smallest imaginary part.
m = 3;
type = 'SM'; % or 'LR' or 'SR', or 'LI', or 'SI'
[V,D] = eigs(A,m,type);
[V,D] = eigs(@(x) MyProd(x,d,v,w),n,m,type);
% For symmetric problems, you can also ask for the m most positive or the
% m most negative. See below.

% 5. It is possible to modify the default parameters that define how eigs
% behaves. New parameter values are set using struct. For example,
% here is how to turn off intermediate printing (?), ease termination
% criteria, reduce the number of iterations, and specify a starting
% vector...

x0 = [1;zeros(n-1,1)];
opts = struct('disp',0,'tol',0.0001,'maxit',50,'v0',x0);
m = 3;
[V,D] = eigs(A,m,'LM',opts);
[V,D] = eigs(@(x) MyProd(x,d,v,w),n,m,'LM',opts);

% 6. To go after non-extremal eigenvalues, you need to supply a function
% that carries out vector = (A - alfa*I)\vector where alfa is a scalar.
% Note that extremal eigenvalues of inv(A - alfa*I) are associated with
% the eigenvalues that are near to alfa. In particular, if tau is an
% eigenvalue of inv(A - alfa*I) then alfa + 1/tau is an eig of A. If
% tau is very big, then this eigenvalue is close to alfa.

n = 20; d = (n:-1:1)'; v = n*randn(n-1,1); w = n*randn(n-1,1);
m = 3;
opts = struct('disp',0);
alfa = 10.1;
[V,D] = eigs(@(x) MyInvProd(x,d,v,w,alfa),n,m,alfa,opts);

% 7. A third output argument can be use to signal lack of convergence
% of one of the eigenvalues...

opts = struct('disp',0,'maxit',3);
d = rand(1000,1); v = randn(999,1); w = randn(999,1);
A = sparseSetUp(d,v,w);
m = 3;
[V,D,flag] = eigs(A,m,'LM',opts);
% if flag==0 then one or more reported eigenvalue untrustworthy.
% Look at residuals norm(A*V(:,k) - D(k,k)*V(:,k)) for more insight.

% 8. If the A-matrix is symmetric, then eigs should be told...

n = 20; d = (n:-1:1)'; v = n*randn(n-1,1);
A = sparseSetUp(d,v,v);
m = 3;
opts = struct('disp',0,'tol',0.0001,'issym',1);
[V,D] = eigs(A,m,'LA',opts);
[V,D] = eigs(@(x) MyProd(x,d,v,v),n,m,'LA',opts);

function A = sparseSetUp(d,v,w)
% d is a column n-vector, v is a column (n-1)-vector, and w is a column
% (n-1)-vector.
% Sets up the n-by-n arrow matrix defined by diagonal d, right column v,
% and bottom row w.
n = length(d);
% A terrible method that sets up the full matrix first...
% A = diag(d); A(1:n-1,n) = v; A(n,1:n-1)=w'; A = sparse(A);

% An element-by-element method...
A = sparse(n,n);
for i=1:n-1
A(i,i) = d(i); A(i,n) = v(i); A(n,i) = w(i);
end
A(n,n) = d(n);

% A vectorized method...
A = sparse(n,n);
i = [1:n 1:n-1 n*ones(1,n-1) ];
j = [1:n n*ones(1,n-1) 1:n-1 ];
s = [ d' v' w' ];
A = sparse(i,j,s);

function y = MyProd(x,d,v,w)
% x is a column n-vector.
% y = A*x where A is the n-by-n arrow matrix defined by diagonal d,
% right column v, and bottom row w.
n = length(x);
y = d.*x + x(n)*[v;0];
y(n) = y(n) + w'*x(1:n-1);

function y = MyInvProd(x,d,v,w,alpha)
% x is a column n-vector.
% y = (A-alpha*I)\x where A is the n-by-n arrow matrix defined by diagonal d,
% right column v, and bottom row w, and alpha is a scalar. Assumes
% that d-alpha is nonzero
n = length(x);
% Look at LU of A - alpha*I. Let D1 = diag(d(1)-alpha,...,d(n-1)-alpha)
% z' = w'*inv(D1), mu = d(n)-alpha...
%
% D1 v I 0 D1 v
% = LU = *
% w' mu z' 1 0 mu-z'*v

mu = d(n)-alpha; d1 = d(1:n-1)-alpha; z = w./d1;
% Overwrite x with L\x..
x(n) = x(n)-z'*x(1:n-1);
% Overwrite x with U\x
x(n) = x(n)/(mu-z'*v); x(1:n-1) = (x(1:n-1)-x(n)*v)./d1;
y = x;