function [X,K,Report] = getXK(DOMAIN,JOB,A,B,G,E,U,V,W,L,RES,Sx,Sr)
% Computes Riccati solution X and gain matrix K from implicit data.
%    DOMAIN       'c' or 'd'.
%    JOB          's' or 'a' (stabilizing/anti-stabilizing)
%    A,B,G,E      Problem data.
%    U,V,W,L,RES  Implicit solver outputs. [U;V;W] spans the (anti)stable
%                 invariant subspace.
%    Sx,Sr        State and R scalings.
%
% On output, REPORT is set to 
%     0    Accurate solution found (success)
%     1    Solution accuracy may be poor
%     2    Solution is not finite
%     3    Pencil has eigenvalues on stability boundary.
%     4    Pencil is singular ([B;S;R] rank deficient)

%   Copyright 2018-2020 The MathWorks, Inc.
Report = 0;  % assume success
if any(isnan(L))
   % [B;S;R] rank deficient
   X = [];  K = [];  Report = 4;
elseif isempty(U)
   X = [];
   K = W;  % watch for static case n=0
   if ~isempty(A)
      Report = 3;
   end
else
   T = matlab.internal.math.nowarn.mrdivide(V,U);  % accuracy limits [eps,1/eps]
   K = -matlab.internal.math.nowarn.mrdivide(W,U);
   nrmT = norm(T,'fro');
   maxRES = max(RES);
   
   if isnan(nrmT) || 10*eps*nrmT>1
      % T is not finite
      Report = 2;
   elseif maxRES>1e-2
      % Residual is too large
      Report = 3;
   else
      % Check closed-loop stability. Solution may not be (anti)stabilizing
      % when
      %   * U nearly singular
      %   * [U;V;W] inaccurate due to eigenvalues near the stability boundary
      % NOTE: When [U;W] is rank deficient, e.g., due to (A,B) not being
      % stabilizable, K=-W/U involves 0/0 divisions and is inaccurate. The
      % computed K may look reasonable but fail to stabilize A-B*K.
      if isempty(E)
         p = eig(A+G*T-B*K);
      else
         p = eig(A+G*T-B*K,E);
      end
      % Check that A-B*K is (anti)stable
      switch [DOMAIN JOB]
         case 'cs'
            STABLE = all(real(p)<0);
         case 'ca'
            STABLE = all(real(p)>0);
         case 'ds'
            STABLE = all(abs(p)<1);
         case 'da'
            STABLE = all(abs(p)>1);
      end
      if ~STABLE
         if sqrt(eps)*nrmT>1
            % U singular is likely cause of failure
            Report = 2;
         else
            % Closed-loop poles near stability boundary are likely cause of failure
            Report = 3;
         end
      end
   end
            
   if Report==0
      % Found finite + stabilizing solution
      if isempty(E)
         X = T;
      else
         X = matlab.internal.math.nowarn.mrdivide(V,E*U);
      end
      % Undo scaling and enforce symmetry
      X = Sx .* X .* Sx';
      K = Sr .* K .* Sx';
      X = (X+X')/2;
      if maxRES>sqrt(eps)
         Report = 1;
      end
   else
      X = [];  K = [];
   end
end