function [s]=schurc(a,x)

%   [s]=schur(a,x)
%   Given an nxn matrix "a" and an order k index vector "x", schurc
%   returns the Schur complement of "a" wrt the principal
%   submatrix indexed by "x". Note that the entries of "x"
%   are automatically put in increasing order.
%         Michael Tsatsomeros 19/11/97

[l,n]=size(a); if l~=n 'matrix is not square',break,end


% Find the complement "y' of "x" in {1,2,...n}, and vectorize 
%              "x" and "y" in increasing order
      y=[]; % initialize y for a vacuous Schur complement
      w=zeros(1,n);
      for k=1:length(x)
           w(x(k))=x(k); %expand x to an 1xn. This automatically
                         %puts entries of "w" in increasing order!
      end  
      z=(1:n)-w;         %create the complement
      i=0;
      for k=1:n          %extract nonzero entries of "z"-->"y"
          if z(k)~=0
             i=i+1; y(i)=z(k);
          end
      end
      i=0;
      for k=1:n          %extract nonzero entries of "w"-->"x"
          if w(k)~=0
          i=i+1; x(i)=w(k);
          end
      end

      c=a(x,x); 
      if det(c)==0 'Schur complement does not exist',break,end
      'The Schur complement wrt', x, 
      s=a(y,y)-a(y,x)*inv(c)*a(x,y);