% define diagonal matrix:
n=1000;

% Case 1
diagA = [1:n]';

% Case 2
%diagA(11:100) = [80*ones(30,1),-90*ones(30,1),-100*ones(30,1)];

% Case 3
% diagA(95:100)=10*diagA(95:100);

A=spdiags(diagA,0,n,n);

% take exact solution and compute righ-hand side b:
xex=ones(n,1);
b=A*xex;

x0 = zeros(n,1); % initial guess

clf % clear figure

for i=1:n/10,
    % make i steps of GMRES:
    [x,H]=gmres_H(A,b,x0,i);

    % plot the residual norm:
    subplot(1,2,1);
    semilogy(i,norm(b-A*x),'*'); hold on
    
    % compute Harmonic Ritz values:
    hrv = eig( (H(1:i,:))'\(H'*H) );

    % plot Harmonic Ritz values and eigenvalues of A:
    subplot(1,2,2);
    plot(hrv,zeros(i,1),'x',diagA,zeros(n,1),'+');
    hold on

    % convert polynomial roots (Harmonic Ritz values)
    % to the GMRES residual polynomial:
    p = poly(hrv);
    p = p/p(i+1);   % take care that p(0)=1

    % plot residual polynomial:
    m = 400; h = (max(diagA)-min(diagA))/m;
    xx = [min(diagA):h:max(diagA)];
    [pp]=polyval(p,xx);
    plot(xx,pp,'-'), hold off

    fprintf('press <Enter>...'); pause, fprintf('\n');
end