%
%  This script illustrates convergence of IRAM with Shift-Invert 
%
%
%  Recommended:  Invoke 
%
%  figure(1); figure(2);
%
%  and resize to taste prior to executing this script
%
%  D.C. Sorensen
%  14 Mar 2000
%  
 format short e
 n = 500;
 
 A = randn(n); f = randn(n,1);
 
 m = 20; k = 5;

%
%  Do m steps of Arnoldi to get an idea of the location of the spectrum of A
%

 H = zeros(m); V = zeros(n,m); [V,R,f] = Arnold(A,V,H,f,1,m);

%
%  Compute the eigenvalues of the projected matrix H
%

 tr = eig(R); k1 = length(tr);
 [s,ii] = sort(-real(tr)); tr = tr(ii);     
 
%
%  Uncomment this statment if you want to see eigenvalues near the 
%
%  RIGHT END of the spectrum
%

 sigma = real(tr(1)) + .1*abs(real(tr(1)) - real(tr(k1))); 
%---------------------------------------------------------

%
%  Uncomment this statment if you want to see eigenvalues near the
%
%  CENTER of the spectrum
%

%sigma = sum(real(tr))/k1;
%--------------------------

 SIGMA = sigma

 f = randn(n,1);


%
%  Compute A1 = (A - sigma I)^{-1}
%
   A1 = inv(A - sigma*eye(n));

%
%  This should NOT BE DONE IN PRACTICE !!!!
%
%  Instead do [L,U] = lu(A - sigma I)  ONCE
%
%  and solve  w = U\(L\v)  when a matrix-vector product 
%
%  is called for in Arnoldi
%


 which1 = 'LM';    
 figure(1)
 [V,R,ritz] = Iram_Plt(A1,which1,k,m,f); 

 R = V'*(A*V);

 Resid__Ritz = [norm(A*V - V*R) norm(ritz)]
 tr = eig(R); k1 = length(tr);
 t  = eig(A);

%
%   sort eigenvalues of A and R by order of distance from the
%   shift sigma
%
    [s,ii] = sort(abs(t-sigma)); t = t(ii);     

    [s,ii] = sort(abs(tr-sigma)); tr = tr(ii);     

%
%        Compute points on a circle centered at sigma 
%        and just enclosing the converged Ritz values
%
         rho = abs(sigma - tr(k1))
         kpts = 50; 
         tau = 2*pi*sqrt(-1)/kpts;
         s = zeros(kpts+1,1);
         for j = 0:kpts,
             s(j+1) = sigma + rho*exp(j*tau);
         end


   eigA__eigR__diff = [ t(1:k1) tr abs(abs(t(1:k1)) - abs(tr))]

%
%  Display eigenvalues of A  with black +'s
%
%  and Ritz values of R  with red o's
%
%  Draw a cyan dashed circle centered at sigma and enclosing
%  the computed Ritz values.
%

 figure(2)
 plot(real(t),imag(t),'+',real(tr),imag(tr),'ro',real(sigma),imag(sigma),'g*',real(s),imag(s),'c--')