function varargout = eigsymqr(A, tol, maxiter)
%symeig Executes the symmetric QR iteration to approximate
%the eigenvalues of a real symmetric matrix A.
%tol is the error tolerance, and maxiter is the maximum number
%of iterations to be performed.
%  E = symeigqr(A) - assigns E a column vector containing
%  the eigenvalues.
%The error is measured by assuring that as the tridiagonal matrix
%T is being reduced to a diagonal matrix,abs(T(k,k-1)) < tol
%in the current k x k submatrix. If the desired tolerance is not
%obtained for the k x k submatrix, the iteration stops and an error
%message is printed.
%Input: symmetric matrix A, error tolerance, tol, and maxiter, the  
%maximum number of iterations to execute.
%Output: column vector E containing the computed eigenvalues.

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

if issym(A) == false
   error('The matrix is not symmetric.');
end

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

[P, T] = trireduce(A);
if nargout == 2
	Q = eye(n);
end

for k=n:-1:2
   I = eye(k);
	iter = 0;
   while abs(T(k,k-1)) > tol*(abs(T(k,k)+ abs(T(k-1,k-1))))
      iter = iter + 1;
		if iter > maxiter
			fprintf('Failure of convergence. ');
			fprintf('Current eigenvalue approximation %g\n',T(k,k));
			break;
		end
		%compute the Wilkinson shift
      r = (T(k-1,k-1)-T(k,k))/2;
		if r ==0
			s = 1;
		else
			s = sign(r);
		end
      sigma = T(k,k)+r-s*sqrt(r^2+T(k,k-1)^2);
      [Q1, R1] = givenshessqr(T(1:k,1:k) - sigma*I);
      T(1:k,1:k) = R1*Q1 + sigma*I;
		if nargout == 2
			P1 = eye(n);
			P1(1:k,1:k) = Q1;
			Q = Q*P1;
		end
   end
   T(k,k-1) = 0;
end

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


