  function [B, Q1] = qrsd(A, nit)
% function  [B  Q]  =  qrsd(A, nit)  
% classical QR  algorithm for computing  the Schur
% form with classical  2x2 shift and deflation the
% only complication of  the code comes from saving
% the Q  matrix.  A smaller  matrix multiplies the
% whole unitary  matrix Q1 - initially  set to Id.
% Q1 is  returned. The  code also ensures  that on
% return  B  =  Q1'  *  A  * Q1  and  B  is  upper
% triangular possibly complex.  Complex arithmetic
% is not avoided.
% 
n = size(A,1) ;
B = A;
Q1 = eye(n,n);
tol = 1.e-12;
%%-------------------- the 'eigenvalue' loop
for ni= n:-1:2
%%-------------------- nit iters. per eigenvalue 
    for i=1:nit
      D=[B(ni-1,ni-1),B(ni-1,ni);B(ni,ni-1),B(ni,ni)];
%%-------------------- shift2 computes the shift
      t =shift2(D);
      fprintf(1,' ni = %d ;   shift = %e\n',ni,t); 
      B(1:ni,1:ni) = B(1:ni,1:ni) - t * eye(ni) ;
      [Q, R]       = qr(B(1:ni,1:ni), 0);
      B(1:ni,1:ni) = R*Q + t*eye(ni)
      B(1:ni,ni+1:n) =   Q'*B(1:ni,ni+1:n);
      B
      pause 
%%-------------------- save transformation Q1
      Q1(:,1:ni)   = Q1(:,1:ni)*Q;
%%-------------------- simple stopping test
      if (norm( B(ni,1:ni-1) ) < tol)
          disp (' *** deflating') 
          break
        end
    end
end
