function varargout = eigb(A,tol,maxiter,recover)
%EIGB Executes the Francis iteration (implicit QR iteration)
%to approximate the eigenvalues and, optionally, the eigenvectors
%of a real matrix A.
%
%   [V D] = eigb(A,tol) returns an orthogonal
%   matrix V and diagonal matrix D of eigenvalues such that
%   V'*A*V = D. tol defaults to 1.0e-8, where tol is the error tolerance.
%
%   E = eigb(A,tol) - assigns E a column vector containing
%   the eigenvalues. The default value of tol is the
%   same as the previous calling sequence.

%   If the desired tolerance is not obtained for any eigenvalue,
%   a warning message is printed and computation continues for the
%   remaining eigenvalues and, optionally, eigenvectors.

[m, n] = size(A);
if m ~= n
	error('The matrix is not square.');
end

if nargin == 1
   tol = 1.0e-8;
   maxiter = 1000;
   recover = 1;
elseif nargin == 2
   maxiter = 1000;
   recover = 1;
elseif nargin ==3
   recover = 1;
end

[P, H] = hess(A);
SH = H;
E = zeros(n,1);
count = 1;
k = n;

if recover
   adjust = 0;
end

iter = 0;

while k > 2
   iter = iter + 1;
   H(1:k,1:k) = impdsqr(H(1:k,1:k));
   if abs(H(k-1,k-2)) < tol
      % there is a complex conjugate pair or two real eigenvalues
      H(k-1,k-2) = 0;
      e = eig(H(k-1:k,k-1:k));
      E(count) = e(1);
      E(count+1) = e(2);
      count = count + 2;
      k = k - 2;
      iter = 0;
   elseif abs(H(k,k-1)) < tol
      % there is one real eigevalue
      H(k,k-1) = 0;
      E(count) = H(k,k);
      count = count + 1;
      k = k - 1;
      iter = 0;
   elseif iter > maxiter
      if recover
         if adjust > 9
            error('The iteration was not successful in computing all the eigenvalues');
         end
         % stir things up and keep trying convergence to lambda_k by doing
         % a single shift iteration
         H(1:k,1:k) = chase(H(1:k,1:k));
         adjust = adjust + 1;
         iter = 0;
      else
         error('The iteration was not successful in computing all the eigenvalues');
      end
   end
end

% compute the final eigenvalues
if k == 2
   e = eig(H(1:2,1:2));
   E(count) = e(1);
   E(count+1) = e(2);
else
   E(count) = H(1,1);
end

if nargout == 2
	%preallocate the n x n matrix of eigenvectors
	V = zeros(n, n);
	%compute the eigenvector for each eigenvalue using the Hessenberg
	%inverse iteration. Apply a maximum of 10 iterations.
	for i=1:n
      x0 = randn(n,1);
      x0 = x0/norm(x0);
	  v = eigvechess(SH, E(i), x0, tol, 10);
      hasinf = isinf(v);
      hasnan = isnan(v);
      if (any(hasinf) == 1) || (any(hasnan) == 1)
         sigma = E(i) + (2+rand)*eps;
         x0 = randn(n,1);
         x0 = x0/norm(x0);
         v = eigvechess(SH,sigma,x0,tol,10);
      end
		V(:,i) = P*v;
	end
end

if nargout == 2
   varargout{1} = V;
   varargout{2} = diag(E);
else
    varargout{1} = E;
end
