function d = qralg(A)
% QR algorithm with shifts and deflation.

% Iterate over eigenvalues
for n = length(A):-1:2
  % QR iteration
  while sum( abs(A(n,1:n-1)) ) > eps
    s = A(n,n);
    [Q,R] = qr(A-s*eye(n));
    A = R*Q + s*eye(n);
  end 
  % Deflation
  d(n) = A(n,n);
  A = A(1:n-1,1:n-1);
end
% Last remaining eigenvalue
d(1) = A(1,1);