% Perform Lanczos with complete reorthogonalization
%
% Inputs:
%   A = input matrix (may be sparse)
%   e = true eigenvalues, sorted from largest to smallest
%   q = starting vector for Lanczos
%   m = number of Lanczos steps to take
%   mm1,mm2 = indices of extreme eigenvalues for which to plot error and error bounds
%   mine, maxe = range in which to plot eigenvalues
%   steplabel = array of values at which to draw vertical lines in plots
%   recomp = 1 to perform algorithm from scratch
%          = 0 just to redisplay data, with possibly different
%               mine, maxe, mm1, mm2, steplabel
%   whichbound = 2 to plot all errors and error bound
%              = 1 to plot local error and error bounds
%              = 0 to plot error bounds only
%
% Outputs:
%   plot of Ritz values for all m Lanczos steps, in range [mine,maxe]
%
%   plot of global error (distance from final eigenvalue) (if whichbound=2)
%           local error (distance from closest eigenvalue) (if whichbound=1)
%           error bound (provided by Lanczos algorithm)
%     for mm1 to mm2 smallest and mm1 to mm2 largest eigenvalues
%
%   plot of components of each eigenvector in current Lanczos vector
%           (only valid if A is diagonal with eigenvalues sorted
%            from largest to smallest!)
%     for mm1 to mm2 smallest and mm1 to mm2 largest eigenvalues
%
if (recomp==1)  % {
   n=max(size(A));
   Eval = zeros(m,m);
   ErrorBounds = zeros(m,m);
   qq = q/norm(q);
   Q = qq;
   clear alpha
   clear beta

   disp('start performing Lanczos, print step number')
   for i=1:m
      if (rem(i,10)==0)
         i,
      end
      z = A*qq;
      alpha(i) = qq'*z;
%     reorthogonalize twice to guarantee orthogonality of Lanczos vectors
      z = z - Q*(Q'*z);
      z = z - Q*(Q'*z);
      beta(i) = norm(z);
      qq = z/beta(i);
      Q = [Q,qq];
   end

   disp('start extracting Ritz values, print step number')
   T = diag(alpha(1:m)) + diag(beta(1:m-1),1) + diag(beta(1:m-1),-1);
   for i=1:m
     if (rem(i,10)==0)
        i,
     end
     [V,D] = eig(T(1:i,1:i));
%    sort eigenvalues into decreasing order
     [Ds,Is] = sort(-diag(D));
     Eval(1:i,i) = -Ds;
     Vs = V(:,Is);
     ErrorBounds(1:i,i) = abs((beta(i)*Vs(i,:))');
   end
   disp('largest deviation of singular value of Lanczos vector matrix Q from 1 ='),
      max(abs(svd(Q)-1))
   disp('done with extracting eigenvalues, bounds')
end % }
%
% Produce graph of converging Ritz values 
%
figure(1)
hold off, clg
for j=1:ceil(i/2),
    rr = 2*j-1;
    if (rem(j,4)==1),
       str='k+';
    elseif (rem(j,4)==2),
       str='r+';
    elseif (rem(j,4)==3),
       str='g+';
    elseif (rem(j,4)==0),
       str='b+';
    end
    deval = diag(Eval,j-1);deval=deval(j:m-j+1);
    plot(rr:m,Eval(j,rr:m),str); hold on
    plot(rr:m,deval,str);
end
plot((m+1)*ones(size(e)),e,'k+')
axis([0,m+2,mine,maxe])
for b=steplabel, plot([b,b],[mine,maxe],'k:'); end
title([int2str(m),' steps of Lanczos (full reorthogonalization) applied to A'])
xlabel('Lanczos step'), ylabel('Eigenvalues')
disp('done with figure 1')
disp('pause (hit return to continue)'),pause
%
% Produce graph of errors, error bounds in largest mm1 to mm2 Ritz values
%
figure(2)
hold off, clg
step = 3;
for js=mm1:step+1:mm2
  hold off, clg
  for j=js:min(js+step,mm2)
    if (rem(j,4)==1),
       str='w';
       strB='--w';
       strD=':w';
    elseif (rem(j,4)==2),
       str='r';
       strB='--r';
       strD=':r';
    elseif (rem(j,4)==3),
       str='g';
       strB='--g';
       strD=':g';
    elseif (rem(j,4)==0),
       str='b';
       strB='--b';
       strD=':b';
    end
    disp(['Eigenvalue being plotted = ',num2str(e(j))]),
%   plot global error
    if (whichbound > 1),
       semilogy((j:m),abs(Eval(j,j:m)-e(j))/(abs(e(j))),str);
       hold on
    end
%   plot local error
    if (whichbound > 0),
       semilogy((j:m), ...
       min(abs(ones(size(e))*Eval(j,j:m)-e*ones(1,m-j+1)))/(abs(e(j))),strD);
       hold on
    end
%   plot error bound
    semilogy((j:m),ErrorBounds(j,j:m)/(abs(e(j))),strB);
    hold on
  end
  axis([0,m+1,1e-16,1])
  jsend = min(js+step,mm2);
  if (whichbound > 1),
   if (jsend ~= js),
      title(['True errors and error bounds, in eigenvalues ',int2str(js), ...
          ' to ', int2str(min(js+step,mm2))])
   else
      title(['True errors and error bounds, in eigenvalue ',int2str(js)])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
   ylabel('Global Error (solid), Local Error (dotted), Bound (dashed)')
  elseif (whichbound > 0),
   if (jsend ~= js),
      title(['True local error and error bounds, in eigenvalues ',int2str(js), ...
          ' to ', int2str(min(js+step,mm2))])
   else
      title(['True local error and error bounds, in eigenvalue ',int2str(js)])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
   ylabel('Local Error (dotted), Bound (dashed)')
  else
   if (jsend ~= js),
      title(['Error bounds, in eigenvalues ',int2str(js), ...
          ' to ', int2str(min(js+step,mm2))])
   else
      title(['Error bounds, in eigenvalue ',int2str(js)])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
  end
  for b=steplabel, plot([b,b],[1e-16,1],'w:'); end
  disp('pause (hit return to continue)'), pause
end
disp('done with figure 2')
disp('pause (hit return to continue)'), pause
%
% Produce graph of errors, error bounds in smallest mm1 to mm2 Ritz values
%
figure(3)
hold off, clg
for js=mm1:step+1:mm2
  hold off, clg
  for j=js:min(js+step,mm2)
    if (rem(j,4)==1),
       str='w';
       strB='--w';
       strD=':w';
    elseif (rem(j,4)==2),
       str='r';
       strB='--r';
       strD=':r';
    elseif (rem(j,4)==3),
       str='g';
       strB='--g';
       strD=':g';
    elseif (rem(j,4)==0),
       str='b';
       strB='--b';
       strD=':b';
    end
    deval = diag(Eval,j-1); deval=deval(1:m-j+1);
    deval = deval';
    disp(['Eigenvalue being plotted = ',num2str(e(n+1-j))]),
%   plot global error
    if (whichbound > 1),
       semilogy((j:m),abs(deval-e(n+1-j))/max(abs(e)),str);
       hold on
    end
%   plot local error
    if (whichbound > 0),
       semilogy((j:m), ...
         min(abs(ones(size(e))*deval-e*ones(1,m-j+1)))/(abs(e(j))),strD);
       hold on
    end
    devalbnd = diag(ErrorBounds,j-1); devalbnd = devalbnd(1:m-j+1);
%   plot error bound
    semilogy((j:m),devalbnd/max(abs(e)),strB);
    hold on
  end
  for b=steplabel, plot([b,b],[1e-16,1],'w:'); end
  axis([0,m+1,1e-16,1])
  jsdif = (n+1-js)-(n+1-min(js+step,mm2));
  if (whichbound > 1),
   if (jsdif ~= 0)
      title(['True errors and error bounds, in eigenvalues ', ...
          int2str(n+1-min(js+step,mm2)),' to ', int2str(n+1-js)])
   else
      title(['True errors and error bounds, in eigenvalue ',  ...
          int2str(n+1-min(js+step,mm2))])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
   ylabel('Global Error (solid), Local Error (dotted), Bound (dashed)')
  elseif (whichbound > 0),
   if (jsdif ~= 0)
      title(['True local error and error bounds, in eigenvalues ', ...
          int2str(n+1-min(js+step,mm2)), ' to ', int2str(n+1-js)])
   else
      title(['True local error and error bounds, in eigenvalue ', ...
          int2str(n+1-min(js+step,mm2))])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
   ylabel('Local Error (dotted), Bound (dashed)')
  else
   if (jsdif ~= 0)
      title(['Error bounds, in eigenvalues ',int2str(n+1-min(js+step,mm2)), ...
          ' to ', int2str(n+1-js)])
   else
      title(['Error bounds, in eigenvalue ',int2str(n+1-min(js+step,mm2))])
   end
   xlabel('Lanczos step (full reorthogonalization)'),
  end
  disp('pause (hit return to continue)'), pause
end
disp('done with figure 3')
disp('pause (hit return to continue)'), pause
%
% Plot components of Lanczos vectors for largest mm1 to mm2 eigenvectors
%
figure(4)
hold off, clg
for js=mm1:step+1:mm2
  hold off, clg
  for j=js:min(js+step,mm2)
    if (rem(j,4)==1),
       str='w';
       strB='--w';
       strQp='+w';
       strQn='ow';
    elseif (rem(j,4)==2),
       str='r';
       strB='--r';
       strQp='+r';
       strQn='or';
    elseif (rem(j,4)==3),
       str='g';
       strB='--g';
       strQp='+g';
       strQn='og';
    elseif (rem(j,4)==0),
       str='b';
       strB='--b';
       strQp='+b';
       strQn='ob';
    end
    disp(['Eigenvalue being plotted = ',num2str(e(j))]),
    Qpos = find(Q(j,1:m)>0);
    Qneg = find(Q(j,1:m)<0);
    semilogy(Qpos,abs(Q(j,Qpos)),strQp);
    hold on
    semilogy(Qneg,abs(Q(j,Qneg)),strQn);
    axis([0,m+1,1e-16,1])
  end
  axis([0,m+1,1e-16,1])
  if (js ~= min(js+step,mm2))
     title(['Lanczos vector components for eigenvectors  ',int2str(js), ...
         ' to ', int2str(min(js+step,mm2))])
  else
     title(['Lanczos vector components for eigenvector  ',int2str(js)])
  end
  xlabel('Lanczos step (full reorthogonalization)'), ylabel('eigencomponent (+ = pos, o = neg)')
  for b=steplabel, plot([b,b],[1e-16,1],'w:'); end
  for b = [1e-12 1e-8 1e-4], plot([0,m+1],[b,b],'w:'); end
  disp('pause (hit return to continue)'), pause
end
disp('done with figure 4')
disp('pause (hit return to continue)'), pause
%
% Plot components of Lanczos vectors for smallest mm1 to mm2 eigenvectors
%
figure(5)
hold off, clg
for js=mm1:step+1:mm2
  hold off, clg
  for j=js:min(js+step,mm2)
    if (rem(j,4)==1),
       str='w';
       strB='--w';
       strQp='+w';
       strQn='ow';
    elseif (rem(j,4)==2),
       str='r';
       strB='--r';
       strQp='+r';
       strQn='or';
    elseif (rem(j,4)==3),
       str='g';
       strB='--g';
       strQp='+g';
       strQn='og';
    elseif (rem(j,4)==0),
       str='b';
       strB='--b';
       strQp='+b';
       strQn='ob';
    end
    disp(['Eigenvalue being plotted = ',num2str(e(n+1-j))]),
    Qpos = find(Q(n+1-j,1:m)>0);
    Qneg = find(Q(n+1-j,1:m)<0);
    semilogy(Qpos,abs(Q(n+1-j,Qpos)),strQp);
    hold on
    semilogy(Qneg,abs(Q(n+1-j,Qneg)),strQn);
    end
  axis([0,m+1,1e-16,1])
  if (n+1-min(js+step,mm2) ~= n+1-js)
     title(['Lanczos vector componentes for eigenvectors ', ...
         int2str(n+1-min(js+step,mm2)), ' to ', int2str(n+1-js)])
  else
     title(['Lanczos vector componentes for eigenvector ',int2str(n+1-min(js+step,mm2))])
  end
  xlabel('Lanczos step (full reorthogonalization)'), ylabel('eigencomponent (+ = pos, o = neg)')
  for b=steplabel, plot([b,b],[1e-16,1],'w:'); end
  for b = [1e-12  1e-8  1e-4 ], plot([0,m+1],[b,b],'w:'); end
  disp('pause (hit return to continue)'), pause
end
disp('done with figure 5')
