function [X,K,L,INFO] = idare(A,B,Q,varargin)
%IDARE  Implicit solver for discrete-time Riccati equations.
%
%   [X,K,L] = IDARE(A,B,Q,R,S,E) computes the stabilizing solution X of
%   the discrete-time algebraic Riccati equation
%                                        -1 
%      E'XE = A'XA - (A'XB + S)(B'XB + R)  (A'XB + S)' + Q
%
%   The matrices Q,R must be Hermitian and E must be nonsingular. When 
%   omitted or set to [], R,S,E default to the values R=I, S=0, E=I. 
%   Scalar-valued Q,R are interpreted as multiples of the identity matrix. 
%   IDARE also returns the state-feedback gain K and the closed-loop
%   eigenvalues L given by
%                    -1
%      K = (B'XB + R)  (B'XA + S'),    L = EIG(A-B*K,E) .
%
%   IDARE returns [] for X,K when there is no finite stabilizing solution.
%
%   [X,K,L,INFO] = IDARE(...) also returns a structure INFO with the
%   following data:
%      * Vectors SX and SR of state and R scalings
%      * Basis [U;V;W] of the stable invariant subspace of the associated
%        scaled matrix pencil
%      * Diagnostic "Report" with value
%           0    Accurate solution found (success)
%           1    Solution accuracy may be poor
%           2    Solution is not finite
%           3    Symplectic spectrum [L;1./L] has eigenvalues on unit circle.
%           4    Pencil is singular ([B;S;R] rank deficient)
%   The solution X and gain K are related to DX=diag(SX), DR=diag(SR), and
%   U,V,W by
%                -1    -1                -1
%      X = DX V U  DX E ,    K = - DR W U  DX .
%
%   [...] = IDARE(...,'noscaling') turns off the built-in scaling and
%   sets all entries of SX and SR to one. This speeds up computation but 
%   can be detrimental to accuracy when A,B,Q,R,S,E are poorly scaled.
%
%   [...] = IDARE(...,'anti') computes the anti-stabilizing solution X
%   that puts all eigenvalues of (A-B*K,E) outside the unit disk.
%
%   See the reference page for additional algorithmic details.
%
%   See also ICARE.

%   Author(s): Pascal Gahinet and Alan J. Laub.
%   Copyright 1986-2018 The MathWorks, Inc.

%   (A-zE,B) must be stabilizable, E invertible, and [B;S;R] full column 
%   rank for a finite stabilizing solution to exist. While these conditions 
%   are not sufficient in general, they become sufficient when R invertible, 
%   [Q S;S' R]>=0, and (A-B*(R\S'),Q-S*(R\S')) detectable.
%
%   Reference: W.F. Arnold, III and A.J. Laub, "Generalized Eigenproblem
%	 Algorithms and Software for Algebraic Riccati Equations," Proc. IEEE, 
%   72(1984), 1746--1754.

% Parse inputs
try
   [A,B,~,Q,R,S,E,AutoScaling,JOB] = ricpack.parseInputs('d',A,B,Q,varargin{:});
catch ME
   throw(ME)
end

% Scale
[n,m] = size(B);
if AutoScaling
   [Sx,Sr,A,B,~,Q,R,S,E] = ricpack.scaleData('d',A,B,[],Q,R,S,E);
else
   Sx = ones(n,1);  Sr = ones(m,1);
end

% Solve in scaled coordinates
[U,V,W,L,RES] = ricpack.DARE(JOB,A,B,Q,R,S,E);

% Compute X,K and construct INFO
[X,K,Report] = ricpack.getXK('d',JOB,A,B,0,E,U,V,W,L,RES,Sx,Sr);
if nargout>3
   INFO = struct('Sx',Sx,'Sr',Sr,'U',U,'V',V,'W',W,'Report',Report);
elseif Report==1
   % Warn if accuracy is poor
   warning(message('Control:foundation:ARE04'))
end