% This script constructs a Krylov matrix K = [b, Ab,..., A^(k-1)b]
% with columns normalized to have length one as in the
% power method sequence.  Then  [V,R] = qr(K,0) is computed
% followed by  H = V'*A*V.
% 
% H should be upper Hessenberg.  This computation illustrates
% that this is not so numerically.

%
%  set up matrix
%
 n = 100; k = 10;

 [Q,R] = qr(randn(n));
 t = max(abs(diag(R)));
 R(n,n) = t*n;
 A = Q*R*Q';
%
%
 v =  randn(n,1); v = v/norm(v);
 K = zeros(n,k);
 K(:,1) = v; 
 resid = zeros(k,1);
%
% Form the Krylov matrix from the power sequence
%
 for j = 1:k,
    w = A*v;
    v = w/norm(w);
    K(:,j) = v;
 end
 format long e
 [V,R] = qr(K,0);
%
% compute the projected matrix H and th residual F
%
 H = V'*(A*V);
 F = A*V - V*H;
 format short e

 Hval = H

 for j = 1:k,
    resid(j) = norm(F(:,j));
 end
%
% plot norms of residual columns
%
 semilogy(resid,'*')
 hold
 semilogy((resid))
 title('Norm of column residuals F = A*V - V*H')
 hold

%
% contour and mesh plots of log abs H
%
 rev = k:-1:1;
 figure
 contour(log(abs(H(rev,:))))    
 colorbar
 title('contour plot log(abs(H))')

 figure
 mesh(log(abs(H(rev,:))))  
 colorbar             
 title('mesh plot log(abs(H))')

%
%  Hcut will show how "Hessenberg" H is
%
 Hcut = H;
 Hnorm = norm(H);
 for j = 1:k,
     for i = 1:k,
         if (abs(Hcut(i,j)) < (10^(-10))*Hnorm), Hcut(i,j) = 0; end
     end
 end
 figure
 spy(Hcut)
 title('Significant part of H should be Hessenberg')