function [V,T,f] = lanczosfplot(A,m,reorthog)
%LANCZOSF computes the Lanczos decomposition of matrix A
%
%  Input:   A -- an n by n symmetric matrix
%           m -- dimension of the Krylov subspace, m << n.
%           reorthog -- 1 if full reorthogonalization and 0 if not
%
%   Output: V -- an n by k orthogonal matrix
%           T -- an m x m symmetric tridiagonal matrix
%           f -- an n vector
%
%           AV = VT + f*e_m', where e_m is the standard basis
%           vector [0 0 ...0 1]'
%           plot of the eigenvalues of T for k = 1, 2, ..., m


%   original author: D.C. Sorensen
%   modified by William Ford
%   15 May 2013

n = size(A,1);
v = rand(n,1);
v1 = v/norm(v);
V = zeros(n,m);
T = zeros(m,m);
eigenvalues = zeros(m,m); 

f = A*v1;
alpha = v1'*f;
f = f - v1*alpha;

V(:,1) = v1;
T(1,1) = alpha;

for j = 2:m,
   beta = norm(f); 
   v0 = v1;
   v1 = f/beta;

   % full reorthogonalization. slow but effective.
   if reorthog == 1
      for i = 1:j-1
         h = V(:,i)'*v1;
         v1 = v1 - V(:,i)*h;
      end
   end

   f = A*v1 - v0*beta;

   alpha = v1'*f;
   f = f - v1*alpha;

   T(j,j-1) = beta;
   T(j-1,j) = beta;
   T(j,j) = alpha;
   V(:,j) = v1;

   % find eigenvalues and eigenvectors for T.
   [~, DT] = eig(T(1:j,1:j));
   sigma = diag(DT);

   % reorder the eigenvalues and the corresponding eigenvectors.
   [~,ind] = sort(sigma,'descend');
   sigma = sigma(ind);   
   for i = 1:j
      eigenvalues(j,i) = sigma(i);
   end
end

plot(eigenvalues, 'o');
