% Antti Hannukainen
%
% Transform A->diag(T+,T-)

% Schur

function [V,D] = tht1(A)

[U,T] = schur(A,'complex');
[US,TS]=ordschur(U,T,'rhp');


% Solve sylvester equation ( triangular system )

% first extract needed matrices :

I = find(real(diag(TS))<0);

n = min(I);
p = size(TS,1);

T11 = TS(1:(n-1),1:(n-1));
T12 = TS(1:(n-1),n:p);
T22 = TS(n:p,n:p);


% Solve Sylvester equation

Z = zeros(size(T11,2),size(T22,2));

for i=1:size(Z,2)
  Z(:,i) = (T11 - T22(i,i)*eye(size(T11)))\( -T12(:,i)+Z*T22(:,i))
end
size(T11)
size(Z)
size(T22)
p-n
n

V = US*[eye(size(T11)) Z ; zeros(p-n+1,n-1) eye(size(T22))];

D = [ T11 zeros(n-1,p-n+1) ; zeros(p-n+1,n-1) T22];
