
% A demonstration of QR iteration, with preliminary reduction,
% to find all eigenvalues and vectors.  H. Motteler, Feb 1999.
%
% This example is intended for transparency and to illustrate
% underlying principles, rather than for efficiency. 

% generate a random matrix A0 w/ real eigenvalues
% (this simple search may fail for n greater than 7 or 8)
n = 7
A0 = rand(n, n);
while ~isreal(eig(A0))
  A0 = rand(n, n);
end

% generate a random symmetric matrix A0
% n = 10
% A0 = rand(n, n);
% A0 = A0'*A0;

A0 = [1 4; 1 1]
n = 2;

% Phase 1: apply Householder similarity transforms to reduce general
% A to Hessenberg form (upper-triangular with one sub-diagonal band),
% or a symmetric A to tridiagonal form.  In either case we will have
% A0 = T'*A*T at the end of this phase.

A = A0;
T = eye(n);
I = eye(n);
for i = 1 : n-2

  % v zeros elements one below the diagonal
  e = zeros(n,1);
  e(i+1) = 1;
  a = A(:,i);
  a(1:i) = zeros(i, 1);
  v = a + sign(a(i+1)) * norm(a) * e ;

  % explicit calculation of H
  H = I - 2 * (v*v') / (v'*v) ;

  % explicit application of similarity transform
  A =  H * A * H ;

  % final T is the product of the transforms H
  T = H * T;

end

% check the preliminary factorization
% norm(T'*A*T - A0)


% Phase 2: do a basic QR iteration with shifts, but without
% deflation.  Note that without deflation, convergence can
% be slow for Hessenberg (initially, general) matrices, and
% the system will occasionally not converge at all.

k = 1;
qr_eps = 1e-4;

% tmask has all ones below diagonal; zero elsewhere
tmask = ((1:n)' * ones(1,n)) > (ones(n,1) * (1:n));

while max(max(abs(tmask .* A))) > qr_eps & k < 1000

  % basic QR iterateion
  % [Q, R] = qr(A);
  % A = R * Q;

  % QR iteration with shifts
  s = A(n,n);
  [Q,R] = qr(A - s*I); 
  A = R * Q + s*I ;

  T = Q' * T;
  k = k + 1;
end

% check final factorization
% norm( T'*A*T - A0 )

% show iteration count
k

% show eigenvalues
% diag(A) 

% compare with Matlab eigenvector routine

norm(sort(diag(A)) - sort(eig(A0)))
