%% Power and inverse power iteration examples
clear all
A=diag([-3 -1 4 5 6],0)
[m,n]=size(A);
T=rand(n);
A = A + T/n % Diagonally dominant
%A=T*A*inv(T) % Similarity transform of diagonal matrix
tol= 1e-10; % Tolerance for calculations
eig(A)

%% Largest eigenvalue using power iteration
[l1,v1,k1] = pow_it(A,tol)

%% Smallest eigenvalue using inverse power iteration
[l4,v4,k4] = inv_pow_it(A,tol)

%% Intermediate eigenvalue using shifts

% Plot Gershgorin disks
nphi=256; t=(0:nphi)*2*pi/nphi;
for i = 1:n
    c = A(i,i); % Centre
    r = sum(abs(A(i,:))) - abs(A(i,i)); % Radius
    plot( r*cos(t)+c, r*sin(t) ,'b-');hold on; axis('equal');grid on; % Disk
    plot(c, 0,'ro','MarkerSize',6,'MarkerFaceColor','red');% Eigenvalue
end
% Plot upper bound disk for all eigenvalues
r=norm(A,inf); % Why inf? Recall that inf is max |row| sum
plot(r*cos(t),r*sin(t),'g-','LineWidth',2);
hold off;

%% Try to find different eigenvalues
% Most negative eigenvalue
shift=-20; [l,v,k] = inv_pow_it(A-shift*eye(n,n),tol); l=shift+l % Shift back

%% Eigenvalue closest to 2
shift=2; [l,v,k] = inv_pow_it(A-shift*eye(n,n),tol); l=shift+l % Shift back

%% What if our initial guess is intermediate between two eigenvalues?
shift=5.5; [l,v,k] = inv_pow_it(A-shift*eye(n,n),tol); l=shift+l 

%% Try a (really) large matrix
n=100;
A=rand(n,n); A=A'*A; % Positive eigenvalues
v=eig(A)
[l1,v1,k1] = pow_it(A,tol)
[ln,vn,kn] = inv_pow_it(A,tol)
%%
l1-v(n)
ln-v(1)