% Illustrate the cubic convergence of tridiagonal QR with Wilkinson's shift
% Inputs:
%   n = dimension of random symmetric matrix to generate
%   m = number of QR steps to take
%
% Generate a randon symmetric n by n matrix, and reduce it to tridiagonal form
% a=randn(n,n);a=a+a';t=hess(a);
disp('Initial tridiagonal matrix')
t
n=min(size(t));
disp('pause (hit return to continue)')
pause
%
% Store last subdiagonal and last diagonal entries of QR iterates in NN
NN = [ t(n,n-1), t(n,n) ];
% Perform QR iteration with Wilkinson's shift
for i=1:m
% Compute the shift
  lc=t(n-1:n,n-1:n);elc=eig(lc);
  if ( abs(t(n,n)-elc(1)) < abs(t(n,n)-elc(2)) ),
      shift = elc(1);
  else
      shift = elc(2);
  end
% Perform QR iteration; enforce symmetry explicitly
  [q,r]=qr(t-shift*eye(n));t=r*q+shift*eye(n);
  t = tril(triu(t,-1),1); t = (t+t')/2;
% Update NN
  NN = [NN;[t(n,n-1),t(n,n)]];
end
NN = [NN,NN(:,2)-NN(m+1,2)]
disp('   T(n,n-1)    T(n,n)     T(n,n)-eigenvalue')
NN
disp('pause (hit return to continue)')
pause
disp('Final tridiagonal matrix')
t