function [G,J,rowcol,rowcol_n]=get_eigen_block_ssts_tssd(M,N,hyp)

% GET_EIGEN_BLOCK_SSTS_TSSD Obtain the eigen decomposition of the Hermitian
% matrix, R, for the generating function such that
%
%      [-T  P   I]
%  R = [ P  0   0]
%      [ I  0   0]
%
% where T is the positive-definite Toeplitz-like Block matrix, and dQ is 
% symmetric, not necessary positive-defininite.
%
% [G,J,rowcol] = GET_EIGEN_BLOCK_SSTS_TSSD(M,N,hyp)
%
% Inputs:
%           M       : Time-series Inputs of Gaussian process, size nx1
%           N       : Time-series Inputs of posterior Gaussian process
%           hyp     : hyperparameters of the covariance function [g v]
%
% Outputs:
%           G       = Generates the Generating function (matrix)
%           J       = J-unitary matrix, J=[I 0;0 -I]
%           rowcol  = block-position 
%
% Note that the covariance function used, is the usual convention, of
%       T(Z|a,g,v) = a * [ exp(-0.5*g*[Zi-Zj]^2) + v ]
%       P(Z|a,g,v) = a * [ exp(-0.5*g*[Zi-tj]^2) ]
%
%
% (C) Gaussian Process SSTS Toolbox 1.0
% (C) Copyright 2006-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

[g_row,e_dec,g_toe,g_dqg,g_to3,g_drn,g_drm]=eigen_tool();

rowcol=g_row(M);
rowcol_n=g_row(N);

n=size(M,1);
Nt=size(N,1);
nn=2*n+Nt;
rk2=size(rowcol,2);
rk2_n=size(rowcol_n,2);
rk=rk2+rk2_n;
A=zeros(rk);
B=zeros(nn,rk);
B(1:n,1)=-g_to3(M,hyp,1);
B(1:n,rk2+1)=g_drn(M,hyp,1,N);
B(n+1:Nt+n,1)=g_drm(M,hyp,1,N);
B(n+Nt+1,1)=1;
for k=2:rk2
    p=rowcol(k);
    tmp1=g_to3(M,hyp,p);
    tmp3=g_to3(M,hyp,p-1);
    tmp2=[0;tmp3(1:end-1,:)];
    B(1:n,k)=tmp2-tmp1;
    tmp1=g_drm(M,hyp,p,N);
    tmp3=g_drm(M,hyp,p-1,N);
    tmp2=[0;tmp3(1:end-1,:)];
    tmp3=tmp1-tmp2;
    B(n+1:Nt+n,k)=tmp3;
end
for k=2:rk2_n
    p=rowcol_n(k);
    tmp1=g_drn(M,hyp,p,N);
    tmp3=g_drn(M,hyp,p-1,N);
    tmp2=[0;tmp3(1:end-1,:)];
    B(1:n,rk2+k)=tmp1-tmp2;
end
tmp0=zeros(nn,rk);
for k=1:rk2
    p=rowcol(k);
    A(k,:)=B(p,:);
    tmp0(p,k)=1;
end
for k=1:rk2_n
    p=rowcol_n(k);
    A(rk2+k,:)=B(n+p,:);
    tmp0(n+p,rk2+k)=1;
end

[G,J]=e_dec(A,B,tmp0);