function E = eigqrshift(A,tol,maxiter)
%eigqrshift shifted Hessenberg QR iteration for computing all the
%eigenvalues of a real matrix whose eigenvalues have distnct magnitudes.
%
%   E = eigqr(A,tol,maxiter), where the tol is error tolerance
%   desired, and maxiter is the maximum number of iterations for
%   computing any single eigenvalue. Vector E contains the eigenvalues
%   of A. If the desired tolerance is not obtained for any particular
%   eigenvalue, a warning message is printed and computation continues
%   for the remaining eigenvalues.


n = size(A,1);
H = hess(A);
normF = norm(H,'fro');

if nargin == 1
   tol = 1.0e-8;
   maxiter = 1000;
elseif nargin == 2
   maxiter = 1000;
end
nonconvergence = false;

for k=n:-1:2
   I = eye(k);
   iter = 0;
   abstrace = abs(H(k,k))+abs(H(k-1,k-1));
   if abstrace ~= 0
      conval = tol*abstrace;
   else
      conval = tol*normF;
   end
	while abs(H(k,k-1)) > conval
      iter = iter + 1;
      if iter > maxiter
			nonconvergence = true;
			break;
      end
      sigma = H(k,k);
      if sigma == 0
         sigma = eps;
      end
      [Q1, R1] = givenshessqr(H(1:k,1:k) - sigma*I);
      H(1:k,1:k) = R1*Q1 + sigma*I;
      abstrace = abs(H(k,k))+abs(H(k-1,k-1));
      if abstrace ~= 0
         conval = tol*abstrace;
      else
         conval = tol*normF;
      end
   end
   H(k,k-1) = 0;
end

if nonconvergence
   disp(['There are eigenvalues that did not converge.'...
         'Consider using another method.']);
end

E = diag(H);
