% QR_eig.m

% function [eig_val,Qtot,iflag] = QR_eig(A,tol,ishow)

%

% This MATLAB file demonstrates the use of iterative QR decompositions

% to obtain the eigenvalues of a matrix.

%

% K.J. Beers. 2/27/2001.

% Version 7/8/2002

function [eig_val,Qtot,iflag] = QR_eig(A,tol,ishow)

iflag = -1;  % signifies no convergence

ITER_MAX = 100;  % maximum # of QR iterations

Awork = A;  % copy of A matrix to work on

if(ne(ishow,0))

    disp('Original matrix A = ');

    disp(A)

end

[Q,R] = qr(Awork);  % perform A = Q*R

eig_val_old = diag(R); % estimate of eigenvalue vector

for iter=1:ITER_MAX

    Awork = R*Q;  % get new matrix similar to A

    [Q,R] = qr(Awork);  % perform QR iteration

    eig_val = diag(R);  % new estimate of eigenvalues

    delta_eig = eig_val - eig_val_old;  % change in eigenvalues

    error = dot(delta_eig,delta_eig);

    eig_val_old = eig_val;  % store old values

    if(mod(iter,ishow)==0)

        disp(['iteration # = ', int2str(iter)]);

        disp(['error = ', num2str(error)])

        disp(Awork)   % display A matrix

    end

    if(error<=tol)

        iflag = 0;  % show convergence achieved

        break;

    end

end

% Display final results.

if(ne(ishow,0))

    disp(' ');

    disp(' ');

    disp('Final results of QR calculation : ');

    disp(' ');

    disp('Eigenvalues : ');

    disp(diag(Awork));

    disp(' ');

    disp('MATLAB eig() eigenvalues : ')

    disp(eig(A));

end

return;