function [U,T] = qrsort(U,T,dt)
%
%  [U,T] = QRSORT(U,T,DT);
%
%  The function QRSORT implements sorting of a Schur decomposition
%  A*U = U*S.
%  As input ORSORT takes the target vector DT (which may be fo any length 
%  up to N where N is the size of the square matrices A,U and S), a unitary
%  matrix U and an upper triangular matrix T.
%  On return, the diagonal entries of T are ordered according to DT.
%  If the #elements in DT are less than N, the diagonal elements of 
%  the original T not spesified in DT follws after those spesified in 
%  DT. The returned (U,T) satisfies: 
%     A = U*T*U', p = #elements in DT, T(i,i) = DT(i) for all i <= p.
%
%  Written by: Kjetil Havre 18/9-1996.
%
%
   [p,m] = size(dt);
   if m > p
      dt = dt.'; h = m; m = p; p = h;
   end
   
   k   = size(T,1);
   nts = max(size(dt));
   if nts < k-1
      t   = diag(T);
      for i = 1:nts
         id = find(t==dt(i));
	 if isempty(id) == 0
            t  = rmvecel(t,id(1));
	 end
      end
      dt = [dt; t];
   end
   
   for i=1:k-1
      t = diag(T(i:k,i:k));
      j = qrselect(dt(i),t)+i-1;
      [U,T] = qrswap(U,T,j,i);
   end
end

