function [U,U_inv,A,Aplus,Aminus,Nminus,pik,lambda] = additive2_decompose_schur(R,Q);

% this function calculates U such that
% U^-1 A U = [0 0 0;0 Aminus 0;0 0 Aplus] 
% such that Aminus has all its eigenvalues in the open left half plane 
% and Aplus has all its eigenvalues in the open right half plane, the first zero is the 
% eigenvalue at the origin
% Here S is the variance matrix, R is the rate matrix, and Q is the
% generator and 
% A = [0 -Q S^-1;
%     I R S^-1];


% S = [10 0; 0 100];
% R = [-4 0;0 2];
% Q = [-1 1;1 -1];
% N is the size of all the matrices
N = length(Q);
e = ones(N,1);
I = eye(N,N);
Z = zeros(N,N);
N;
size(R);
%S_inv = S;
%for i=1:N
%    S_inv(i,i) = 1/S(i,i);
%end
R_inv = R;
for i=1:N
    R_inv(i,i) = 1/R(i,i);
end

% pik is the stationary vector for Q
%pik = null(Q')';
%size(pik)
%size(e)
%pik = pik/(pik*e);
Q1=[Q ones(N,1)];
pik=[zeros(1,N) 1]/Q1;
% lambda is the mean drift
lambda = pik*R*e;
if(abs(lambda)<10^-5)
    disp('WARNING!!!-lambda is zero-WARNING!!!')
else
%if(lambda<0)
%    Nminus = N;
%    Nplus = N-1;
%else
%    Nminus = N-1;
%    Nplus = N;
%end
Nminus=-1;
Nplus=0;
for i=1:N
    if R(i,i)>0
        Nminus=Nminus+1;
    else
        Nplus=Nplus+1;
    end
end
if(lambda<0)
    Nminus = Nminus+1;
    Nplus = Nplus-1;
end
A=Q*R_inv;
[Q2,R2] = schur(A); 
z = -10000;
[Q1,R1] = SRSchur(Q2,R2,z,Nminus+1);
    A11 = R1(1:Nminus+1,1:Nminus+1);
    A11;
    A12 = R1(1:Nminus+1,Nminus+2:N);
    A22 = R1(Nminus+2:N,Nminus+2:N);
    A22;
    X1 = lyap(A11,-A22,A12);
    X11 = X1(1:Nminus,:);
    X12 = X1(Nminus+1,:);
    Aminus = A11(1:Nminus,1:Nminus);
    Aminuszero = A11(1:Nminus,Nminus+1);
    X2 = -Aminus\Aminuszero;
    temp = [eye(Nminus) X2 X11;zeros(1,Nminus) 1 X12;zeros(Nplus,Nminus+1) eye(Nplus)];
    temp_inv = [eye(Nminus) -X2 -X11+X2*X12;zeros(1,Nminus) 1 -X12;zeros(Nplus,Nminus+1) eye(Nplus)];
    U = Q1*temp;
    U_inv = temp_inv*Q1';
    % U_inv = inv(U);
    Aplus = A22;
    U;
    U_inv*A*U;
    Aplus;
    Aminus;
end