function [Q,T] = reordschr(Q,T,dt,epp);
%
%  [U,T] = reordschr(U,T,DT,EPP);
%
%  Assumes that (U,T) is a complex schur form of the matrix A = U*T*U';
%  where T is strictly upper triangular. 
%  Then the returned (U,T) satisfies: 
%     A = U*T*U', p = #elements in DT, T(i,i) = DT(i) for all i <= p.
%  EPP - tolerance used to check if T(i,i) == DT(i) 
%
%  Written by: Kjetil Havre 18/9-1996.
%
%  Example sorting eigenvalues and eigenvectors in ascending order:
%
%  [U1,T1] = schur(A);
%  [U,T] = rsf2csf(U1,T1);
%  P = diag(T2);
%  Ps = sort(P);
%  [U,T] = reordschr(U,T,Ps);
%
%  Note that QRSORT does the same but much faster. 
%
   if nargin < 4
      epp = 1e-14;
   end
   
   [p,m] = size(dt);
   if m > p
      dt = dt.'; h = m; m = p; p = h;
   end
   
   n = max(size(Q));
   t = diag(T); ts = t(1:p);
   
   ndf = norm(dt - ts,inf);
   if ndf < epp
%      disp('Already in desired form')
      return;
   end
   
   ncdt = max(size(dt));
   if ncdt > 1   
      for i = 1:ncdt
         [Q,T] = reordschr(Q,T,dt(ncdt-i+1),epp);
      end
      return
   end
      
   while ndf > epp,
      for k=1:n-1
          itk  = find( abs(dt-t(k)) < epp );
	  itk1 = find( abs(dt-t(k+1)) < epp );
          if isempty(itk) & (isempty(itk1)==0)
	     G = givens( T(k,k+1), T(k+1,k+1) - T(k,k) );
	     T(k:k+1,k:n) = G*T(k:k+1,k:n);
	     T(1:k+1,k:k+1) = T(1:k+1,k:k+1)*G';
	     Q(1:n,k:k+1) = Q(1:n,k:k+1)*G';
             t = diag(T); 
	     ts = t(1:p);
	  end
          ndf = norm(dt - ts,inf);
      end
   end   
end
