function [B,V,k] = qr_it(A,tol,kmax);
% Usage: [B,V,k] = qr_it(A,tol,kmax);
% Computes all eigenpairs with iteration tolerance tol by
% QR iteration.
% Input:
% A = matrix
% tol = tolerance
% kmax = max allowable number of iterations
% Output:
% l = vector of eigenvalues
% V = matrix of associated eigenvectors
% k = number of iterations
% err = error at each iteration


[n,n]=size(A);
V = eye(n);
k = 1;
err = norm(A);
while err>tol & k<kmax
    A_old = A;
    [Q,R] = qr(A); % Compute QR factorization
    
    A = R*Q;       % Generate next matrix
    V = V*Q;       % Update product of unitary transforms
    err = norm(A-A_old);
    k = k+1;
end
if k>kmax
    error('Not converged')
end
B=A;