% Matlab code qriter.m
% For "Applied Numerical Linear Algebra",  Figures 4.3, 4.4
% Written by James Demmel, Oct 22, 1995
%                Modified, Jun  5, 1997
%
% Illustrate convergence of unshifted QR iteration
% Used to produce figures 4.2 and 4.3
%
% Inputs:
%  D = vector of eigenvalues 
%  m = number of iterations
%
% Outputs:
%  plots of convergence of Orthogonal Iteration (or QR iteration)
%  applied to a matrix with eigenvalues D and randomly chosen
%  eigenvectors
%
[Ds,Is] = sort(-abs(D));
Ds = D(Is);
n = length(D);
Ds = reshape(Ds,n,1);
S = randn(n,n);
disp(['condition number(eigenvector matrix) = ',num2str(cond(S))])
A = S*diag(Ds)*inv(S);

Asave=A;
econv=0; clear econv;
rconv=0; clear rconv;
econv(1:n,1) = abs( diag(A) - Ds );
for j=1:n
  rconv(j,1) = norm(A(j+1:n,1:j));
end
disp('Starting matrix='), Asave
disp('(pause - hit return to continue)'), pause
for i=1:m+1
  [Q,R]=qr(A);
  Alast = A;
  A = R*Q;
  disp(['Matrix ',int2str(i),' = ']),A
  disp('(pause - hit return to continue)'), pause
  econv(1:n,i) = abs( diag(A) + sort(-Ds) );
  for j=1:n-1,
    rconv(j,i) = norm(A(j+1:n,1:j));
  end
end
figure(1), hold off, clf
figure(2), hold off, clf
sbplt = 0;
for i = 1:n,
   if ( rem(sbplt,4) == 0 ) & ( i ~= 1 ),
      disp('hit return to continue'), pause,
      figure(1), hold off, clf
      figure(2), hold off, clf
      sbplt = 0;
   end
   sbplt = sbplt+1;
   figure(1)
   subplot(2,2,sbplt)
   if i == 1,
     rat = abs(Ds(2)/Ds(1));
   elseif i == n,
     rat = abs(Ds(n)/Ds(n-1));
   else
     rat = max(abs(Ds(i)/Ds(i-1)),abs(Ds(i+1)/Ds(i)));
   end
   dconv = (rat*ones(m+1,1)).^((0:m)');
   semilogy(econv(i,:),'-b'), grid, hold on
   semilogy(dconv,'--b'), 
   title(['Convergence of lambda(',int2str(i),')= ',num2str(Ds(i))])
   xlabel('actual error (solid), estimate (dashed)')
   if ( i < n )
     figure(2)
     subplot(2,2,sbplt)
     rrat = abs(Ds(i+1)/Ds(i));
%    rconv = (rrat*ones(m+1,1)).^((0:m)');
%    semilogy(econv(i,:),'-b'), grid, hold on
%    semilogy(dconv,'--b'), 
     dconv = (rrat*ones(m+1,1)).^((0:m)');
     semilogy(rconv(i,:),'-b'), grid, hold on
     semilogy(dconv,'--b'), 
     title(['Convergence of Schur form(',int2str(i+1),':',int2str(n),' , ', ...
                                         '1:',int2str(i),')']),
     xlabel('actual error (solid), estimate (dashed)')
   end
end