function [invQ,logdet,tr,invB]=schur_invssts(G,ord,b)

% SCHUR_INVSSTS       Generates the Generalised Schur's complement of
% the Hermitian matrix, R, such that
%
% R=[-T  I]
%   [ I  0]
%
% where T is the positive-definite Toeplitz-like Block matrix. Note
% that no hyperbolic rotation is implemented.
%
% [Q,logdet,tr,invB] = SCHUR_INVSSTS(G,ord,b)
%
% Inputs:
%           G   : generating matrix, size (2n)x(2rk)
%           ord : order of matrix
%           b   : vector for inv(T)*b
%
% Outputs:
%           Q       = returns N*inv(T)*N'
%           logdet  = log(det(T))
%           tr      = trace(inv(T))
%           invB    = returns column vector of inv(T)*b
% 
% such that N is a permuting matrix, defined by ord.
%
%
% (C) Gaussian Process SSTS Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

n=size(G,1);
n2=n/2;
np=n2+1;
na=np+1;
rk=size(G,2)/2;
d0=ceil(n2/(ord+1));
ind=(ord+1)*ones(d0,1);
ind(1,1)=1;
ind=cumsum(ind);
dx=1;
invQ=zeros(d0);
invB=zeros(n2,1);
logdet=0;
tr=0;
for k=1:n2
    nk=na-k;
    g=G(1,:);
    l=G*[g(1,1:rk)'; -g(1,rk+1:end)'];
    d=l(1,1);
    logdet=logdet+log(abs(d));
    
    if dx<d0
        if k==ind(dx+1)
            dx=dx+1;
        end
    end
    u=l(nk:np,1);
    for r=1:dx
        tx=u(ind(r))/d;
        for s=1:r
            invQ(r,s)=invQ(r,s)-tx*u(ind(s));
        end
    end
 
    r=(l(nk:np,1)'*b(1:k))/d;
    invB(1:k)=invB(1:k)-l(nk:np,1)*r;
    r=l(nk:np,1)'*l(nk:np,1);
    tr=tr-r/d;
    
    tmp=l*(g/d);
    G=G+blaschke_potapov(tmp,n2);
    G=G(2:end,:);      
end
for r=2:d0
    for s=1:(r-1)
        invQ(s,r)=invQ(r,s);
    end
end


function X=blaschke_potapov(Gmp,n2)
N=size(Gmp,1);
rk=size(Gmp,2);

tmp=[zeros(1,rk);Gmp(1:end-1,:)];
tmp(N-n2+1,:)=zeros(1,rk);
X=tmp-Gmp;