% Matlab code tridiQR.m
% For "Applied Numerical Linear Algebra",  Chapter 5, Example 5.7
% Written by James Demmel, Nov 13, 1995
%                 Modified Jun  6, 1997
%
% Illustrate the cubic convergence of tridiagonal QR with Wilkinson's shift
% Inputs:
%   t = symmetric tridiagonal matrix on which to run QR
%   m = number of QR steps to take
%
% You can also generate a random symmetric n by n matrix, 
% and reduce it to tridiagonal form as follows:
%
%    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('   ')
disp('    T(n,n-1)    T(n,n)       T(n,n)-eigenvalue')
disp('    --------    ------       -----------------')
NN
disp('Final tridiagonal matrix')
t