## Copyright (C) 1996, 1998 Auburn University.  All rights reserved.
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## Octave is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
## for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{Xinf},@var{x_ha_err}] =} hinfsyn_ric(@var{A},@var{BB}.@var{C1},@var{d1dot},@var{R},@var{ptol})
## @table
## @item forms
##        xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1 BB'];
##        Ha = [A 0*A; -C1'*C1 -A'] - xx;
## @end table
## and solves associated Riccati equation
## returns error code
## @table
## @item x_ha_err
##    0: successful
##    1: Xinf has imaginary eigenvalues
##    2: Hx not Hamiltonian
##    3: Xinf has inf. eigenvalues (numerical overflow)
##    4: Xinf not symmetric
##    5: Xinf not positive definite
##    6: R is singular
## @end table
## @end deftypefn

function [Xinf, x_ha_err] = hinfsyn_ric (A, BB, C1, d1dot, R, ptol)

  x_ha_err = 0;        # assume success
  Xinf = [];                 # default return value
  n = is_square(A);
  nw = is_square(R);
  if(rank(R) != nw)    x_ha_err = 6;
  else                 # build hamiltonian Ha for X_inf
    xx = ([BB; -C1'*d1dot]/R) * [d1dot'*C1, BB'];
    Ha = [A, 0*A; -C1'*C1, -A'] - xx;
    x_ha_err = 0;
    [d, Ha] = balance(Ha);
    [u, s] = schur(Ha, "A");
    rev = real(eig(s));

    if (any(abs(rev) <= ptol))  # eigenvalues near the imaginary axis
      x_ha_err = 1;
    elseif (sum(rev > 0) != sum(rev < 0))
      ## unequal number of positive and negative eigenvalues
      x_ha_err = 2;
    else
      ## compute positive Riccati equation solution
      u = d * u;
      Xinf = u(n+1:2*n,1:n) / u(1:n,1:n);
      if (!all(all(finite(Xinf))))
        x_ha_err = 3;
      elseif (norm(Xinf-Xinf') >= 10*ptol)
        ## solution not symmetric
        x_ha_err = 4;
      else
        ## positive semidefinite?
        ## force symmetry (faster, avoids some convergence problems)
        Xinf = (Xinf + Xinf')/2;
        rev = eig(Xinf);
        if (any(rev <= -ptol))
          x_ha_err = 5;
        endif
      endif
    endif
  endif
endfunction