function [V,D,nit] = QRA_sym_diag(A,tol)
% function [D,V] = QRA_sym_diag(A,acc)
% QR algorithm without shift for diagonalizing a symmetric matrix
% tol - tolerance for numerical zero, by default, tol = max(size(A)) * eps(norm(A))

if nargin==1
    tol = max(size(A)) * eps(norm(A));
end
if norm(A-A')>tol
    error('Matrix is not symmetric!');
end
[V,D]=hess(A);
nit=0;
while norm(D-diag(diag(D)))>tol
    [Q,R]=qr(D);
    D=R*Q;
    V=V*Q;
    nit=nit+1;
end
D=diag(diag(D));
end

