%IDREIG_DRIVER is an example driver for the IDREig package.
%
%   The script IDREig_driver generates some matrices and
%   vectors and calls five different functions:
%
%     1) IDRORes_academic, an academic implementation of
%   the prototype IDR(s) [Sonneveld/van Gijzen, SISC 31, 2008].
%     2) IDRORes, the less academic implementation of the same
%   prototype IDR(s), which generates only the basic output
%   of IDRORes_academic, but saves workspace.
%     3) IDREig, a routine that presents how to extract
%   eigenvalue information from a run of IDRORes(_academic)
%   based on five different approaches outlined in the report
%   [Gutknecht/Zemke, 2010].
%     4) IDRRitzPurified, a routine that shows how to cheaply
%   estimate the accuracy of the eigenpair approximation.
%     5) compute_W, a routine that exemplifies how stable the
%   computation of a certain set of basis vectors useful for
%   the computation of certain Ritz vectors is for the chosen
%   matrix and starting vector.
%
%   See also: IDRORes, IDRORes_academic, IDREig, IDRRitzPurified.

%   Coypright 2009-2010 Jens-Peter M. Zemke.
%   Version 1.00, Date: 2010/04/13 14:20:34 CEST.

n = 50;
A = randn(n)+(1+sqrt(n))*eye(n);
b = A*ones(n,1);
x0 = randn(n,1);
s = 3;
PH = orth(randn(n,s)+i*randn(n,s)).';
maxj = 20;
start = 1;

[x,resvec,uH,UT,C,omega,trueresvec,R] = ...
                    IDRORes_academic(A,b,x0,s,PH,maxj,start);
[x2,resvec2,uH2,UT2,C2,omega2] = ...
                    IDRORes(A,b,x0,s,PH,maxj,start);

if sum(abs([norm(x-x2) norm(resvec-resvec2) ...
            norm(full(uH-uH2)) norm(full(UT-UT2)) ...
            norm(C-C2) norm(omega-omega2)])) ~= 0
  error('Something went terribly wrong.')
end

colordef black;

figure(1)
steps = (s+1)*maxj;
semilogy(0:steps,trueresvec,'b-',...
         0:steps,resvec,'r-.');
axen = axis;
axis([0 steps axen(3:4)])
legend('true residual norm',...
       'computed residual norm')

[uSH,SR,uSM,PR,uDH,DR,uLM] = IDREig(uH,UT,C,omega,s,maxj);

SH = full(uSH(1:size(SR,1),:));
SM = full(uSM(1:size(SR,1),:));
DH = full(uDH(1:size(DR,1),:));
LM = full(uLM(1:size(DR,1),:));

if issparse(A)
  eigA = eigs(A,min(200,n/2),'SM');
else
  eigA = eig(A);
end

eigSP = eig(SH,full(SR));
eigSM = eig(SM);
eigPP = eig(SH,full(PR));
eigDP = eig(DH,full(DR));
eigLM = eig(LM);

figure(2)
plot(real(eigA),imag(eigA),'r+',...
     real(eigSP),imag(eigSP),'mo',...
     real(eigSM),imag(eigSM),'ms',...
     real(eigPP),imag(eigPP),'gd',...
     real(eigDP),imag(eigDP),'bv',...
     real(eigLM),imag(eigLM),'b^',...
     real(1./omega),imag(1./omega),'wx',...
     real(eigA),imag(eigA),'r+');
legend('eigenvalues of A',...
       'eigenvalues of the Sonneveld pencil',...
       'eigenvalues of the Sonneveld matrix',...
       'eigenvalues of the purified pencil',...
       'eigenvalues of the deflated pencil',...
       'eigenvalues of the deflated matrix',...
       'inverse local minimizers')

if maxj*(s+1) < 301

  lambda = max(eigA);
  [v,theta,acc,estconv,trueconv,allconv,backerr] = ...
     IDRRitzPurified(A,R,uSH,SR,PR,omega,resvec,lambda,s,maxj);
  
  normA = sqrt(norm(A,inf)*norm(A,1));
  figure(3)
  semilogy(1:steps,trueconv/normA,'b-',...
           1:steps,estconv/normA,'r-',...
           1:steps,backerr/normA,'y-',...
           1:steps,allconv/normA,'g-',...
           1:steps,trueconv/normA,'b-',...
           1:steps,estconv/normA,'r-',...
           1:steps,backerr/normA,'y-');
  axis([1 steps 1e-16 1e5])
  legend('true distance to eigenvalue',...
         'estimated rate of convergence',...
         'backward error of the Ritz pair',...
         'distances to other Ritz values')
  title(sprintf('convergence to the eigenvalue lambda = %g',lambda));

end

[RW,pnRW,W,pnW] = compute_W(A,LM,R(:,1),omega,resvec,s,maxj);

figure(4)
semilogy(1:steps,resvec(1:steps),'r-.',...
         1:steps,trueresvec(1:steps),'b-',...
         1:steps,pnW,'w-',...
         1:steps,pnRW,'g-',...
         1:steps,abs(resvec(1:steps)-pnRW'),'r.')
legend('updated residual norms',...
       'true residual norms',...
       'purified residual norms',...
       'recomputed residual norms',...
       'differences')
top = 10^(round(max(log10(resvec))));
axis([1,steps,1e-21*top,1e2*top]);
