% Illustrate convergence of unshifted QR iteration
% Used to produce figures 4.2 and 4.3
% Copyright 1995 by James Demmel
%
% Inputs:
%  D = vector of eigenvalues 
%  m = number of iterations
%
% Outputs:
%  plots of convergence of Orthogonal Iteration, or QR iteration
%
[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
for i=1:m+1
  [Q,R]=qr(A);
  A = R*Q;
  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, clg
figure(2), hold off, clg
sbplt = 0;
for i = 1:n,
   if ( rem(sbplt,4) == 0 ) & ( i ~= 1 ),
      disp('hit return to continue'), pause,
      figure(1), hold off, clg
      figure(2), hold off, clg
      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,:),'-w'), grid, hold on
   semilogy(dconv,'--w'), 
   title(['Convergence of lambda(',int2str(i),')= ',num2str(Ds(i))])
   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,:),'-w'), grid, hold on
     semilogy(dconv,'--w'), 
     title(['Convergence of Schur form(',int2str(i+1),':',int2str(n),' , ', ...
                                         '1:',int2str(i),')']),
   end
end