## Copyright (C) 1996 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301 USA.

## -*- texinfo -*-
## @deftypefn {Function File} {} zgreduce (@var{sys}, @var{meps})
## Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren,
## Automatica, # 1982).
## @end deftypefn

function retsys = zgreduce (Asys, meps)

  ## SYS_INTERNAL accesses members of system data structure

  is_digital(Asys);             # make sure it's pure digital/continuous

  exit_1 = 0;                   # exit_1 = 1 or 2 on exit of loop

  if(Asys.n + Asys.nz == 0)
    exit_1 = 2;                 # there are no finite zeros
  endif

  while (! exit_1)
    [Q,R,Pi] = qr(Asys.d);              # compress rows of D
    Asys.d = Q'*Asys.d;
    Asys.c = Q'*Asys.c;

    ## check row norms of Asys.d
    [sig,tau] = zgrownorm(Asys.d,meps);

    ## disp("=======================================")
    ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ...
    ##   ", tau=",num2str(tau)])
    ## sysout(Asys)

    if(tau == 0)
      exit_1 = 1;               # exit_1 - reduction complete and correct
    else
      Cb = Db = [];
      if(sig)
        Cb = Asys.c(1:sig,:);
        Db = Asys.d(1:sig,:);
      endif
      Ct =Asys.c(sig+(1:tau),:);

      ## compress columns of Ct
      [pp,nn] = size(Ct);
      rvec = nn:-1:1;
      [V,Sj,Pi] = qr(Ct');
      V = V(:,rvec);
      [rho,gnu] = zgrownorm(Sj,meps);

      ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)])
      ## Cb
      ## Db
      ## Ct
      ## Sj'

      if(rho == 0)
        exit_1 = 1;     # exit_1 - reduction complete and correct
      elseif(gnu == 0)
        exit_1 = 2;     # there are no zeros at all
      else
        mu = rho + sig;

        ## update system with Q
        M = [Asys.a , Asys.b ];
        [nn,mm] = size(Asys.b);

        pp = rows(Asys.d);
        Vm =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
        if(sig)
          M = [M; Cb, Db];
          Vs =[V',zeros(nn,sig) ; zeros(sig,nn), eye(sig)];
        else
          Vs = V';
        endif
        ## disp("zgreduce: before transform: M=");
        ## M
        ## Vs
        ## Vm

        M = Vs*M*Vm;

        ## disp("zgreduce: after transform: M=");
        ## M

        ## disp("debugging code:")
        ## Mtmp = [Asys.a Asys.b; Asys.c Asys.d]
        ## Vl = [V', zeros(nn,mm); zeros(mm,nn),Q]
        ## Vr =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)];
        ## Mtmpf = Vl*Mtmp*Vr

        idx = 1:gnu;
        jdx = nn + (1:mm);
        sdx = gnu + (1:mu);

        Asys.a = M(idx,idx);
        Asys.b = M(idx,jdx);
        Asys.c = M(sdx,idx);
        Asys.d = M(sdx,jdx);

        ## disp(["zgreduce: resulting system: nn =",num2str(nn)," mu=",num2str(mu)])
        ## sysout(Asys)
        ## idx
        ## jdx
        ## sdx
      endif
    endif
  endwhile

  ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]);

  if(exit_1 == 2)
    ## there are no zeros at all!
    Asys.a = Asys.b = Asys.c = [];
  endif

  ## update dimensions
  if(is_digital(Asys))
    Asys.nz = rows(Asys.a);
  else
    Asys.n = rows(Asys.a);
  endif

  retsys = Asys;

endfunction