function [A1,k] = my_qr_eig(A,tau,tol,kmax)

% Input:  A : an n by n matrix;
%        tol: tolerance
%        tau: shift factor
%
% Output: A1: An approximate diagonal matrix whose diagonal elements are 
%            eigenvalues of A.
%         k: Number of iterations.
% Note: This code is best used for symmetric matrix.

error = 1; A0 = A;   [m n] = size(A);k=0;
while (error >= tol & k < kmax)
   k = k +1;
   [q r] = qr(A0-tau*eye(n));
   A1 = r*q + tau*eye(n); 
   error = 0;
   for i=3:n;
     for j=1:i-2;
 	error = error + abs(A1(i,j));
     end
   end 
   A0 = A1;
end