function [U,U_inv,A,Aplus,Aminus,Nminus,pik,lambda] = additive3_decompose_ordered_schur(Splus,Rplus,Rstar,Qplusplus,Qpluszero,Qplusstar,Qzeroplus,Qzerozero,Qzerostar,Qstarplus,Qstarzero,Qstarstar,sminus);
% 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
[splus,sstar] = size(Qplusstar);
[splus,szero] = size(Qpluszero);
e=ones(splus+sstar+szero,1);
% pik is the stationary vector for Q
pik = null([Qplusplus Qpluszero Qplusstar;Qzeroplus Qzerozero Qzerostar;Qstarplus Qstarzero Qstarstar]')';
pik = pik/(pik*e);

% lambda is the mean drift
lambda = pik*[Rplus zeros(splus,sstar+szero);zeros(szero,szero+splus+sstar);zeros(sstar,splus+szero) Rstar]*e;
N=sstar+splus+szero;
if(lambda<0)    
    Nminus = splus+(sstar-sminus);
    Nplus = splus+sminus-1;
else   
    Nminus = splus+(sstar-sminus)-1;
    Nplus = splus+sminus;
end
    
    A=[-Rplus zeros(splus,sstar) eye(splus);...
        Splus zeros(splus,sstar) zeros(splus);...
        zeros(sstar,splus) -Rstar zeros(sstar,splus)];
    B=[ (-Qplusplus+Qpluszero*Qzerozero^-1*Qzeroplus) (-Qplusstar+Qpluszero*Qzerozero^-1*Qzerostar) zeros(splus);...
        zeros(splus) zeros(splus,sstar) eye(splus);...       
        (-Qstarplus+Qstarzero*Qzerozero^-1*Qzeroplus) (-Qstarstar+Qstarzero*Qzerozero^-1*Qzerostar) zeros(sstar,splus)];
    [Q2,R2] = schur(B*A^-1);
    
    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:2*splus+sstar);
    A22 = R1(Nminus+2:2*splus+sstar,Nminus+2:2*splus+sstar);
    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;
    