function [e]=ppt(a,x)

%   [e]=ppt(a,x)
%   Given an nxn matrix "a" and an order k index vector "x", 
%   ppt returns the principal pivot transform of "a" relative 
%   to the principal submatrix indexed by "x"
%    (The entries of "x" are put in ascending order)


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

% Find the complement "y" of "x" in <n>, and vectorize 
% "x" and "y" in increasing order
      % initialize y for the possibility of a vacuous Schur complement
      y=[];
      w=zeros(1,n);
      for k=1:length(x)
           w(x(k))=x(k); % expand x to an 1xn. This automatically
                         % puts the 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" into "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" into "x"
          if w(k)~=0
          i=i+1; x(i)=w(k);
          end
      end
c=a(x,x); d=inv(c);
%if det(c)==0 'p.p.t. does not exist',break,end
'The p.p.t. relative to', x, 
e(x,x)=d;
e(x,y)=-d*a(x,y);
e(y,x)=a(y,x)*d;
e(y,y)=a(y,y)+a(y,x)*e(x,y);