function [Q,R,ap] = SRSchur(Q,R,z,b);

% SYNTAX: [Q,R,ap] = SRSchur(Q,R,z,b)
%
% INPUT: orthogonal real Q and quasi-triangular real R such that AQ=QR and a 
% target z in the complex plane. The fourth parameter b determines the length  
% of the ordering with respect to z to be produced:
%
% if b < 0 then -b blocks will be sorted,
% if b > 0 then  b or b+1 eigenvalues will be sorted, depending on the sizes 
% of the blocks,
% if b = 0 then the whole Schur form will be sorted.
%
% OUTPUT: orthogonal real Q and quasi-triangular real R such that AQ=QR with 
% the diagonal blocks ordered with respect to the target z. The number of
% ordered blocks/eigenvalues is determined by the parameter b.
% A vector ap warns for inaccuracy of the solution if an entry of ap exceeds
% one. 
% 
% SUBFUNCTIONS: normalize.m, swaplist.m, select.m, swap.m, lu_complpiv.m
%
% SEE ALSO: schur.m, rsf2csf.m

r = find(abs(diag(R,-1)) > 100*eps);
s = 1:size(R,1)+1;
s(r+1) = [];
%mm=cputime;
for k=1:length(s)-1;
  sk = s(k);
  if s(k+1)-sk == 2
    [Q,R] = normalize(Q,R,sk:s(k+1)-1);
    p(k)  = R(sk,sk)+sqrt(R(sk+1,sk)*R(sk,sk+1));
  else
    p(k)  = R(s(k),s(k));
  end
end
for k = swaplist(p,s,z,b); 
  v      = s(k):s(k+1)-1; 
  w      = s(k+1):s(k+2)-1;
  nrA    = norm(R([v,w],[v,w]),inf);
  [Q,R]  = swap(Q,R,v,w); 
  s(k+1) = s(k)+s(k+2)-s(k+1);
  v      = s(k):s(k+1)-1; 
  w      = s(k+1):s(k+2)-1;
  if length(v)==2
    [Q,R] = normalize(Q,R,v);
  end
  if length(w)==2
    [Q,R] = normalize(Q,R,w);
  end
  ap(k)  = norm(R(w,v),inf)/(10*eps*nrA);
end
%cputime-mm
%disp('3')
R = R - tril(R,-2);
for j=2:length(s)-1; R(s(j),s(j)-1)=0; end
%cputime-mm
%disp('4')
% ----------------------------------------------%

  function [U,S] = normalize(U,S,v);
  n  = size(S,1);
  Q  = rot(S(v,v));
  S(:,v) = S(:,v)*Q;
  S(v,:) = Q'*S(v,:); 
  U(:,v) = U(:,v)*Q;

  % ----------------------------------------------%

  function Q = rot(X);
  c = 1; s = 0;
  if X(1,1)~=X(2,2);
    tau   = (X(1,2)+X(2,1))/(X(1,1)-X(2,2));
    off   = sqrt(tau^2+1);
    v     = [tau - off, tau + off];
    [d,w] = min(abs(v));
    c     = 1/sqrt(1+v(w)^2);
    s     = v(w)*c;
  end
  Q = [c -s;s c];

% ----------------------------------------------%

  function v = swaplist(p,s,z,b);
  n = length(p);
  k = 0; v = [];
  srtd = 0; 
  q = diff(s);
  fini = 0;
  while ~fini
    k        = k+1;
    [dum,j]  = select(p(k:n),z);
    p(k:n+1) = [p(j+k-1) p(k:n)];
    p(j+k)   = [];
    q(k:n+1) = [q(j+k-1) q(k:n)];
    q(j+k)   = [];
    v        = [v,j+k-2:-1:k];
    srtd     = srtd + q(k);
    fini     = (k==n-1)|(k==-b)|(srtd==b)|((srtd==b+1)&(b~=0));
  end

% ----------------------------------------------%

  function [val,pos] = select(p,z);
  y = real(z)+abs(imag(z))*i;
  [val pos] = min(abs(p-y));

% -----------------------------------------------%

  function [U,S] = swap(U,S,v,w);  
  [p,q] = size(S(v,w)); Ip = eye(p); Iq = eye(q);
  r = [];
  for j=1:q
    r = [r;S(v,w(j))];
  end 
  K = kron(Iq,S(v,v))-kron(S(w,w)',Ip);
  [L,H,P,Q] = lu_complpiv(K);
  e = min(abs(diag(H)));
  sigp = 1:p*q;
  for k = 1:p*q-1;
    sigp([k,P(k)]) = sigp([P(k),k]);
  end 
  r = e*r(sigp); 
  x = (H\(L\r));
  sigq = 1:p*q;
  for k = 1:p*q-1;
    sigq([k,Q(k)]) = sigq([Q(k),k]);
  end
  x(sigq) = x;
  X = [];
  for j=1:q
    X = [X,x((j-1)*p+1:j*p)];
  end
  [Q,R]      = qr([-X;e*Iq]);
  S(:,[v,w]) = S(:,[v,w])*Q;
  S([v,w],:) = Q'*S([v,w],:); 
  U(:,[v,w]) = U(:,[v,w])*Q;

  % -----------------------------------------------%

    function [L,U,P,Q] = lu_complpiv(A);
    P = []; Q = []; n = size(A,1);
    for k=1:n-1;
      [a,r] = max(abs(A(k:n,k:n)));
      [dummy,c] = max(abs(a));
      cl  = c+k-1;
      rw  = r(c)+k-1;
      A([k,rw],:) = A([rw,k],:);
      A(:,[k,cl]) = A(:,[cl,k]);
      P(k) = rw; Q(k) = cl;
      if A(k,k) ~= 0;
        rs = k+1:n;
        A(rs,k)  = A(rs,k)/A(k,k);
        A(rs,rs) = A(rs,rs)-A(rs,k)*A(k,rs);
      end
    end
    U = tril(A')'; L = tril(A,-1) + eye(n);
   
  

  
