% This script demonstrates loss of orthogonality in the Arnoldi
 % process (using CGS orthogonalization)
 %
 %  set up matrix
 %
  n = 200; k = 10;
 
  [Q,R] = qr(randn(n));
  t = max(abs(diag(R)));
  R(n,n) = t*100;
  A = Q*R*Q';
 %
 %
  v =  randn(n,1); 
 [V,H,f] = Arnoldi(A,k,v);

 ek = zeros(k,1); ek(k)= 1;
 Resid = norm(A*V - V*H -  f*ek');


 orthtestV = norm(eye(k) - V'*V);
 ortestf = norm(V'*f);

    disp ('     k          Resid      orthf         orthV  ')
    disp([ k    Resid ortestf orthtestV ])

 t = eig(A);
 t = sort(abs(t));

 s = sort(abs(eig(H)));
 abs_eigA_eigH = [t(n-k+1:n) s]


 while (k < n),
    disp('strike a key')
    pause
    k = k + 10 
  
   
    [V,H,f] = Arnoldi(A,k,v);
   
    ek = zeros(k,1); ek(k)= 1;
    Resid = norm(A*V - V*H -  f*ek');
   
   
   
    orthtestV = norm(eye(k) - V'*V);
    ortestf = norm(V'*f);
    disp ('     k          Resid      orthf         orthV  ')
    disp([ k    Resid ortestf orthtestV ])
   
    t = eig(A);
    t = sort(abs(t));
   
    s = sort(abs(eig(H)));
    k1 = 10;
    abs_eigA_eigH__diff = [t(n-k1+1:n) s(k-k1+1:k)  abs(t(n-k1+1:n) - s(k-k1+1:k)) ]

end