%% QR iteration example
clear all;
% Distinct eigenvalues, symmetric matrix
A=[-1 2 3 -2 5; 2 4 -2 3 1; 3 -2 1 5 7; -2 3 5 2 -9; 5 1 7 -9 3]

% Non-distinct eigenvalues
%A=diag([-1 -1 1 3 4 5 -5 5]);[m,n]=size(A);T=rand(n);A=inv(T)*A*T

% Reference solution from matlab
[V_exact,D] = eig(A)
% Sort by magnitude of eigenvalues
l_exact=diag(D);
[y,idx]=sort(abs(l_exact),'descend')
l_exact = l_exact(idx);
V_exact=V_exact(:,idx);

l_exact
V_exact


%% Find all eigenpairs using QR iteration
tol=1e-8; kmax=1000;
t=cputime;[B,V,k] = qr_it(A,tol,kmax);tend=cputime-t;

disp('Rayleigh Quotient iteration');
disp(['Number of iterations = ' sprintf('%d',k) ' , cpu time = ', sprintf('%0.5e',tend)]);
disp(' ');

B 
l=diag(B);
disp('Eigenvalues'); diag(B)
disp('Eigenvectors'); V

disp(['Eigenvalue error = ', sprintf('%0.5e',norm(abs(l)-abs(l_exact))), ' Eigenvector error = ', sprintf('%0.5e',norm(abs(V)-abs(V_exact)))]);

%% Find all eigenpairs using QR iteration with shifts
tol=1e-8; kmax=1000;
t=cputime;[l,V,k] = qr_shift(A,tol,kmax);tend=cputime-t;

% Sort by magnitude of eigenvalues
[y,idx]=sort(abs(l),'descend');
l = l(idx);
V=V(:,idx);


disp('Rayleigh Quotient iteration with shifts');
disp(['Number of iterations = ' sprintf('%d',k) ' , cpu time = ', sprintf('%0.5e',tend)]);
disp(' ');
disp('Eigenvalues'); l
disp('Eigenvectors'); V;
disp(['Eigenvalue error = ', sprintf('%0.5e',norm(abs(l)-abs(l_exact))), ' Eigenvector error = ', sprintf('%0.5e',norm(abs(V)-abs(V_exact)))]);