function [X, L, G] = dare (A, B, Q, R, S=[], E=[])

##
## function [X, L, G] = dare (A, B, Q, R, S, E)
##
## finds stabilizing solution of discrete algebraic Riccati equation
##
## E'XE = Q + A'*X*A - (A'*X*B + S)*inv(R + B'*X*B)*(B'*X*A + S');
##
## Optimal feedback law is: u = -G*x
##
##
## Method: generalized eigenvalue problem as described in:
##
## Pappas, Laub, Sandell, IEEE TAC, AC-25, 4, 1980, 631-641.
## (Handles singular A)
##
## Emami-Naeini, Franklin, IEEE TAC AC-25, 5, 1980, 1015-1016.
## Ionescu, Weiss, Linear Algebra and Its Applications, 174: 229-238  (1992).
## (Handles singular A and singular R)
##
## example: A and R both singular
## A=[0, 1; 0, -1]; B=[1, 0; 2, 1]; Q=[-4, -4; -4, 7]/11; S=[3,1;-1,7]; R=[9, 3; 3, 1];
##   solution: 
## X =[ -1.4021    13.0569; 13.0569  -125.6365]
## L =[-0.21706 0.68727]
## G =[0.94045  -11.00984; -1.78286   19.60900]
##
## Note: algorithm doesn't know/care if [Q, S; S', R] \geq 0
## but guarantees eigenvalues of (A-BG)  are stable
##
##  James B. Rawlings, 12/31/2018
##
## Argument processing adapted from:
## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
## Created: October 2009
## Version: 0.5.1

  if (nargin < 4 || nargin > 6)
    print_usage ();
  endif

  if (! is_real_square_matrix (A, Q, R))
    ## error ("dare: a, q, r must be real and square");
    error ("dare: ##s, ##s, ##s must be real and square", ...
            inputname (1), inputname (3), inputname (4));
  endif
  
  if (! is_real_matrix (B) || rows (A) != rows (B))
    ## error ("dare: a and b must have the same number of rows");
    error ("dare: ##s and ##s must have the same number of rows", ...
            inputname (1), inputname (2));
  endif
  
  if (columns (R) != columns (B))
    ## error ("dare: b and r must have the same number of columns");
    error ("dare: ##s and ##s must have the same number of columns", ...
            inputname (2), inputname (4));
  endif

  if (! is_real_matrix (S) && ! size_equal (S, B))
    ## error ("dare: s(##dx##d) must be real and identically dimensioned with b(##dx##d)",
    ##         rows (s), columns (s), rows (b), columns (b));
    error ("dare: ##s(##dx##d) must be real and identically dimensioned with ##s(##dx##d)", ...
            inputname (5), rows (s), columns (s), inputname (2), rows (B), columns (B));
  endif

  if (! isempty (E) && (! is_real_square_matrix (E) || ! size_equal (E, A)))
    ## error ("dare: a and e must have the same number of rows");
    error ("dare: ##s and ##s must have the same number of rows", ...
            inputname (1), inputname (6));
  endif

  ## check stabilizability
  if (! isstabilizable (A, B, E, [], 1))
    ## error ("dare: (a, b) not stabilizable");
    error ("dare: (##s, ##s) not stabilizable", ...
            inputname (1), inputname (2));
  endif

  [n,m] = size(B);

  if (isempty(S))
    S = zeros(n,m);
  endif
  if (isempty(E))
    E = eye(n);
  endif

## set up generalized eigenvalue problem Mx = lambda Nx
  M = [A, zeros(n), B; ...
       Q, -E, S; ...
       S', zeros(m,n), R];
 
  N = [E, zeros(n), zeros(n,m); ...
       zeros(n), -A', zeros(n,m); ...
       zeros(m,n), -B', zeros(m)];

## ordered Q,Z factorization with stable eigenvalues first
  [MM, NN, Z, LAMBDA] = qz(M, N, "S");

## stable?
  if (max(abs(LAMBDA(1:n))) > 1 )
    warning("dare: stabilizing solution does not exist")
  endif

  Z1 = Z(1:n, 1:n);
  if (rcond(Z1) < 100*eps)
    warning("dare: badly conditioned or stabilizing solution does not exist")
  endif
  Z2 = Z(n+1:end, 1:n);
  XG = Z1' \ Z2';
  X = XG(:,1:n);
  G = -XG(:,n+1:end)';
  L = eig(A - B*G);
endfunction
