 function [Q,T] = subsit0(A, Q, nit) 
%% [Q,T] = subsit(A, Q, nit) 
%% does "nit" steps of subspace iteration 
%% (with a Schur projection every 10 steps)
%% on return
%% Q   = orth. basis of invariant subspace
%% T   = Q'*A*Q = projection of A on Q.
%%------------------------------------------------ 
 Q = qr(Q,0); AQ = A*Q;
 for i=1:nit
   [Q,R] = qr(AQ,0); 
   AQ = A*Q;
   if (mod(i,10) == 0) 
     [U,T] = schur(Q'*AQ);
%%     [U,T] = eig(Q'*AQ);
     Q = Q*U;
     AQ  = A*Q;  
     R   = AQ - Q*T; 
     res =  norm(R,'fro');
    fprintf(1,' it %3d ; res: %8.3e ',i,res);
    fprintf(1,' T11,T22 %10.8f %10.8f\n',T(1,1),T(2,2)) 
   end 
 end


