function [invQ,NQf]=schur_invssts_cz_matlab(G,ord,n)

% SCHUR_INVSSTS_CZ_MATLAB       Generates the fast algorithm for finding
% the standard deviation of the State-Space Time-Series Gaussian Process
% Prior Model using Generalised Schur algorithm. Two outputs are obtained
% from this aspect.
%
%      [-T   A   I]
%  R = [ A'  0   0]
%      [ I   0   0]
%
% where T is positive-definite strongly Toeplitz-like Block and dQ is the
% cross-covariance matrix between the explanatory variables and the
% predicted explanatory variables.
%
% [invQ,NQf] = SCHUR_INVSSTS_CZ_MATLAB(G,ord,n)
%
% Inputs:
%           G   : generating matrix, size (N0+Nt+1)x(rk+rk_n+1)x2
%           ord : order of transformation
%           n   : original datasize
%
% Outputs:
%           invQ    = returns diagonal of A'*inv(T)*A
%           NQf     = returns the product of N*inv(T)*A
%
%
% (C) Gaussian Process SSTS Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

ng=size(G,1);
rk=size(G,2)/2;
np=ng-2*n;
ord1=ord+1;
d0=ceil(n/ord1);
invQ=zeros(np,1);
NQf=zeros(d0,np);
dw=1;
for k=1:n
    g=G(1,:);
    l=G*[g(1,1:rk)'; -g(1,rk+1:end)'];
    d=l(1,1);
    
    dx=n-k+2;
    dy=ng-n-k+1;
    tmp1=l(dx:dy)/d;
    invQ=invQ-l(dx:dy).*tmp1;
    dx=0;
    dy=ng-n-k+2;
    for h=1:dw
        NQf(h,:)=NQf(h,:)-tmp1'*l(dy+dx);
        dx=dx+ord1;
    end
    if k==(ord1*dw)
        dw=dw+1;
    end

    tmp=l*(g/d);
    G=G+blaschke_potapov_3(tmp,n,np);
    G=G(2:end,:);
end


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

% function GG=hyper_rotate(G);
% 
% rk=size(G,2)/2;
% J=[eye(rk) zeros(rk);zeros(rk) -eye(rk)];
% p=G(1,:)*J*G(1,:)';
% clear J;
% GG=zeros(size(G));
% 
% x=G(1,1:rk);
% a=norm(x);
% g1=x+[a zeros(1,rk-1)];
% q1=eye(rk)-2/(g1*g1')*g1'*g1;
% GG(:,1:rk)=G(:,1:rk)*q1;
% 
% x=G(1,rk+1:rk*2);
% a=norm(x);
% g1=x+[zeros(1,rk-1) a];
% q1=eye(rk)-2/(g1*g1')*g1'*g1;
% GG(:,rk+1:rk*2)=G(:,rk+1:rk*2)*q1;
% 
% clear a g1 q1 x;
% 
% sx=[GG(:,1) GG(:,rk*2)];
% sa=sx(1,1); sb=sx(1,2);
% 
% G3=sx*[1 1;-1 1];
% if p>0
%    G3=G3*(0.5*sqrt([(sa+sb)/(sa-sb) 0;0 (sa-sb)/(sa+sb)]));
% else
%    G3=G3*(0.5*sqrt([(sa+sb)/(sb-sa) 0;0 (sb-sa)/(sa+sb)]));
% end
% G3=G3*[1 -1;1 1];
% 
% GG(:,1)=G3(:,1);
% GG(:,2*rk)=G3(:,2);