function varargout = dare(A,B,Q,varargin)
%DARE  Solve discrete-time algebraic Riccati equations.
%
%   DARE is obsolete, please use IDARE instead.
%
%   See also IDARE.

% Old help:
%DARE  Solve discrete-time algebraic Riccati equations.
%
%   [X,L,G] = DARE(A,B,Q,R,S,E) computes the unique 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
%
%   or, equivalently (if R is nonsingular)
%                                -1             -1                 -1
%    E'XE = F'XF - F'XB(B'XB + R)  B'XF + Q - SR  S'  with  F:=A-BR  S'.
%
%   When omitted, R, S and E are set to the default values R=I, S=0,
%   and E=I.  Beside the solution X, DARE also returns the gain matrix
%                        -1
%          G = (B'XB + R)  (B'XA + S'),
%
%   and the vector L of closed-loop eigenvalues (i.e., EIG(A-B*G,E)).
%
%   [X,L,G,REPORT] = DARE(...) returns a diagnosis REPORT with
%   value:
%     * -1 if the Symplectic matrix has eigenvalues on the unit circle
%     * -2 if there is no finite stabilizing solution X
%     * Frobenius norm of relative residual if X exists and is finite.
%   This syntax does not issue any error message when X fails to exist.
%
%   [X1,X2,D,L] = DARE(A,B,Q,...,'factor') returns two matrices 
%   X1, X2 and a diagonal scaling matrix D such that X = D*(X2/X1)*D.
%   The vector L contains the closed-loop eigenvalues.  All outputs 
%   are empty when the associated Symplectic matrix has eigenvalues 
%   on the unit circle.
%
%   See also CARE.

%   Author(s): Alan J. Laub (1993) with key contributions by Pascal Gahinet, 
%              Cleve Moler, and Andy Potvin
%   Revised: 94-10-29, 95-07-20, 95-07-24, 96-01-09
%   Copyright 1986-2008 The MathWorks, Inc.

%       Assumptions: E is nonsingular, Q=Q', R=R', and the associated
%                    symplectic pencil has no eigenvalues on the unit circle.
%       Sufficient conditions to guarantee the above are stabilizability,
%       detectability, and [Q S;S' R] >= 0.

%	Reference: W.F. Arnold, III and A.J. Laub, ``Generalized Eigenproblem
%	           Algorithms and Software for Algebraic Riccati Equations,''
%	           Proc. IEEE, 72(1984), 1746--1754.

% Flags and backward compatibility
varargin = controllib.internal.util.hString2Char(varargin);
isFlag = cellfun('isclass',varargin,'char');
Flags = varargin(:,isFlag);
IMPLICIT = any(strcmp(Flags,'implicit'));
REPORT = any(strcmp(Flags,'report'));
FACTORS = IMPLICIT || any(strcmp(Flags,'factor'));
if any(strcmp(Flags,'nobalance'))
   NewFlags = {'noscaling'};
else
   NewFlags = {};
end

try
   % Process numeric arguments
   [A,B,~,Q,R,S,E] = ricpack.parseInputs('d',A,B,Q,varargin{~isFlag});
   % Call IDARE
   [X,G,L,INFO] = idare(A,B,Q,R,S,E,NewFlags{:});
catch ME
   throw(ME)
end

% Set REPORT flag
switch INFO.Report
   case {3,4}
      Report = -1;  % Hamiltonian pencil has jw-axis eigenvalues
   case 2
      Report = -2;  % there is no finite stabilizing solution X
   otherwise
      Report = 0;
end

% Process results and build outputs
if FACTORS
   % Factored form
   X1 = INFO.U;
   X2 = INFO.V;
   if isempty(X1)
      D = [];
   else
      Sx = INFO.Sx;
      D = diag(Sx);
      nx = size(A,1);
      if IMPLICIT
         if isempty(E)
            Z = orth([Sx .\ X1 ; Sx .* X2]);
         else
            Z = orth([E * (Sx .\ X1) ; Sx .* X2]);
         end
         X1 = Z(1:nx,:);   X2 = Z(nx+1:2*nx,:);
      elseif ~isempty(E)
         Z = orth([Sx .* E * (Sx .\ X1); X2]);
         X1 = Z(1:nx,:);   X2 = Z(nx+1:2*nx,:);
      end
   end
   if IMPLICIT
      varargout = {X1 X2 L Report};
   else
      varargout = {X1 X2 D L Report};
   end
else
   no = nargout;
   if no<4 && ~REPORT
      switch Report
         case -1
            error(message('Control:foundation:ARE22'))
         case -2
            error(message('Control:foundation:ARE20'))
      end
   elseif no==4 && Report==0
      % Return residual in REPORT
      if isempty(E)
         E = 1;
      end
      T1 = A'*X*A - E'*X*E;
      T2 = (A'*X*B + S)*G;
      Res = T1 - T2 + Q;
      Report = norm(Res,1)/(1+norm(T1,1)+norm(T2,1)+norm(Q,1));
   end
   varargout = {X L G Report};
end