function varargout = eigsymb(A, tol, maxiter)
%eigsymb Executes the Francis iteration (implicit QR iteration)
%to approximate the eigenvalues and, optionally, the eigenvectors
%of a real symmetric matrix A.
%
%   [V D] = eigsymb(A,tol,maxiter) returns an orthogonal
%   matrix V and diagonal matrix D of eigenvalues such that
%   V'*A*V = D. tol defaults to 1.0e-12, and maxiter to 50.
%   tol is the error tolerance, and maxiter is the maximum number
%   of iterations to be performed.
%
%   E = symeigqr(A,tol,maxiter) - assigns E a column vector containing
%   the eigenvalues. The default values of tol and maxiter are 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 issym(A) == false
   error('The matrix is not symmetric.');
end

if nargin == 1 || nargin == 0
   tol = 1.0e-10;
   maxiter = 500;
elseif nargin == 2
   maxiter = 500;
end

if nargout == 2
   eigenvectors = true;
else
   eigenvectors = false;
end

[P, T] = trireduce(A);
if eigenvectors == true
   Q = eye(n);
end
for k=n:-1:2
	iter = 0;
   while abs(T(k,k-1)) > tol*(abs(T(k,k)+ abs(T(k-1,k-1))))
		if iter > maxiter
			fprintf('Failure of convergence. ');
			fprintf('Current eigenvalue approximation %g\n',T(k,k));
			break;
      end
      if eigenvectors == true
         [Q1, T(1:k,1:k)] = chase(T(1:k,1:k),n);
         Q = Q*Q1;
      else
         T(1:k,1:k) = chase(T(1:k,1:k),n);
      end
		iter = iter + 1;
   end
   T(k,k-1) = 0;
	T(k-1,k) = 0;
end

T = diag(diag(T));
if nargout == 2
	P = P*Q;
   varargout{1} = P;
   varargout{2} = T;
else
    varargout{1} = diag(T);
end
