   function [Q,T,res] = subsit1(A, Q, sig, nit, p) 
%% function [Q,T,res] = subsit1(A, Q, sig, nit, p) 
%%
%% does "nit" steps of subspace iteration 
%% (with a Schur projection every 10 steps)
%% on return
%% uses B = A - sig*I to iterate    
%% Q   = orth. basis of invariant subspace
%% T   = Q'*A*Q = projection of A on Q.
%%------------------------------------------------ 
%% Q = qr(Q,0); AQ = A*Q;
 k = 0;
 for i=1:nit
   Q = A*Q - sig*Q;
   if (mod(i,10) == 0)
     [Q,~] = qr(Q,0); 
     AQ = A*Q;
     [U,T] = eig(Q'*AQ);
%%-------------------- select p dominant eigenvalues
     [~, indx] = sort(abs(diag(T)),'descend');
     indx = indx(1:p);
%%-------------------- check 'invariance' measure
     R = AQ(:,indx) - Q(:,indx)*T(indx,indx);
     t =  norm(R,1);
     k = k+1;   res(k) = t; 
     fprintf(1,' it %3d ; res: %10.4e  \n',i,t);
%%-------------------- new Q for next iteration
     Q = Q*U;
   end 
 end


